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TECHNICAL  REPORT  SUMMARY 


One  purpose  of  this  investigation  has  been  to  regionalize 
the  upper  several  hundred  kilometers  of  the  mantle  under  the 
Arctic  region,  Siberia  and  the  Eurasian  continental  area  using 
seismic  surface  waves.  A second  purpose  has  been  to  develop  and 
test  efficient  computer  techniques  for  the  computation  of  accurate 
theoretical  seismograms  for  hypothetical  sources  located  within 
Eurasia  and  which  make  use  of  the  results  of  the  regionalization 
part  of  this  study.  These  theoretical  seismograms  are  complete 
in  the  sense  that  they  contain  both  body  and  surface  waves;  they 
can  be  applied  directly  in  a discrimination  program  by  comparing 
the  theoretical  seismograms,  computed  for  both  earthquakes  and 
underground  explosions,  with  the  actual  recorded  seismogram. 

Regionali zation . The  program  of  regionalization  involves 
the  measurement  of  phase  velocities  for  Rayleigh  waves  travers- 
ing the  regions  under  investigation.  The  phase  velocity  curves 
were  obtained  with  the  single-station  phase  velocity  method, 
which  is  described  in  the  main  text  of  this  report.  Briefly, 
the  method  involves  selecting  records  of  earthquakes  from  the 
library  of  the  World  Wide  Standardized  Seismographic  Network 
(WWSSN)  for  events  which:  (1)  occurred  within,  or  on  the  peri- 
meter of,  the  area  of  interest,  (2)  produced  good  long-period 
surface  wave  records  at  WWSSN  stations  for  which  the 
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epicenter-to-station  lines  lie  within  the  regions  being  in- 
vestigated, and  (3)  generated  good,  large  recordings  at  a suf- 
ficient number  of  stations  to  ensure  an  accurate  fault  plane 
solution.  When  a suitable  earthquake  was  found,  an  extensive 
data  processing  and  data  reduction  system  was  applied  to  trans- 
form the  data  into  phase  velocity  curves  for  each  of  the 
selected  epicenter-to-station  lines. 

The  complete  set  of  phase  velocity  curves  (for  fundamental- 

i 

mode  Rayleigh  waves  recorded  by  the  WWSSN  instruments) , which  j 

j 

are  then  used  for  the  regionalization  of  the  Arctic  region,  j 

i 

Siberia,  and  the  Eurasian  continental  area,  consists  of  about  j 

i 

50  dispersion  curves.  Based  on  this  data  set,  we  have  regionalized  i 

I 

the  area  of  interest  by  making  use  of  an  assumption  that  large  j 

parts  of  the  region  with  similar  geophysics  or  basement  geology  j 

j 

and  age  will  have  similar  upper  mantle  structures.  A similar  • 

i 

assumption  was  made  with  success  in  the  regionalization  of  the  j 

i 

Pacific  Ocean  basin  (Kausel,  Leeds  and  Knopoff,  1974;  Leeds,  | 

Knopoff  and  Kausel,  1974). 

The  inversion  has  proceded  by  the  use  of  both  linearized 
and  non-linear  inversions  procedures  with  the  mantle  properties  j 

in  each  of  the  geographic  regions  as  unknowns.  Our  regionalization  I 

has  made  use  of  maps  of  basement  geology  of  Eurasia  coupled  with  j 

f 

a postulate  of  similarity  of  mantle  cross-section  for  regions  of  j 

) 

similar  basement  age,  made  previously  in  small-scale  regional  j 

studies  (Knopoff,  1972) . j 


From  our  regionalization  studies  on  the  Eurasian  continent, 
we  draw  the  following  conclusions;  First,  the  properties  of  the 
mantle  of  the  sediment-covered  Siberian  and  Russian  platforms 
are  highly  consistent  with  those  of  the  Baltic  and  Siberian 
shields  and  with  ancient  shields  observed  worldwide  (Knopoff, 

1972)  . Second,  the  Tibetan  plateau  has  an  extremely  thick 
crust,  perhaps  as  great  as  75  km  from  surface  to  Moho.  Further- 
more, the  upper  mantle  under  the  Tibetan  plateau  has  high  seismic 
velocities,  indicating  the  absence  of  partial  melting  to  relatively 
great  depths.  This  can  be  accounted  for  by  emplacement  of  the 
Indian  shield  under  Tibet  during  the  collision  of  the  Asian  and 
Indian  plates.  Third  the  Alpide  folded  belt  of  Iran  and  Turkey 
has  a very  well-developed  low-velocity  channel  in  the  mantle, 
with  good  contrast  to  the  lid  above,  implying  the  presence  of  a 
zone  of  partial  melting  at  a depth  of  about  90  to  100  km.  Fourth, 
the  Mongolia  -Sinkiang  geophysical  province  has  a well-developed 
low-velocity  zone  in  the  upper  mantle  and  the  more-or-less  stable 
part  of  Eastern  China  also  has  a low  velocity  channel  at  typical 
depths . 

The  principal  problem  in  the  inversion  has  been  the  con- 
struction of  the  boundaries  to  the  geophysical  provinces,  for  which 
only  incomplete  information  is  found  in  the  literature.  A full 
resolution  of  this  problem  is  still  in  the  future,  but  the 
regionalization  used  to  date,  in  which  Eurasia  is  divided  into  6 
rather  large  regions, is  statistical] consistent  with  the  data 
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and  has  about  the  correct  number  of  degrees  of  freedom  in  the 
model  parameterization. 

Theoretical  Seismograms.  The  second  phase  of  our  work  has 
been  the  development  and  testing  of  techniques  for  the  computation 
of  realistic  theoretical  seismograms  at  the  WWSSN  and  HGLP  in- 
stallations,  for  earthquake  and  explosive  sources  within  the  areas 
under  investigation.  For  SH  waves  excited  by  both  earthquakes 
and  explosions,  we  can  construct  complete  seismograms  down  to  a 
period  of  10  sec.  (and  occasionally  to  shorter  periods) . The 
difficulty  of  calculating  diffraction  effects  due  to  the  presence 
of  lateral  heterogeneity  remains  as  an  important  unsolved  problem; 
as  long  as  the  waves  cross  regional  boundaries  at  conditions  re- 
mote from  grazing  incidence,  our  theoretical  seismograms  are 
probably  quite  accurate.  The  extension,  development  and 
optimization  of  these  algorithms  and  computer  programs  for 
Rayleigh,  or  P-SV  waves  on  a spherical,  gravitating  earth  has 
also  been  one  of  the  main  efforts  under  this  contract.  The 
efficient  construction  of  accurate,  theoretical  SH  seismograms  for 

realistic  models  of  the  earth's  structure  using  multimode  methods 
f 

I ■ is  based  on  effective  dispersion  computations,  attenuation 

calculations,  structural  transformations,  and  computations  of 
eigenfunction  characteristics,  effects  of  sphericity,  and 
point-source  response;  we  have  studied  all  these  ingredients  in 


I 


detail. 

Our  work  with  theoretical  seismograms  has  concerned  the 
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multi-mode  surface  wave  phases  Lg,  Sa  and  Sn;  references  to 
this  work  are  appended.  For  discrimination  studies,  it  is  also 
desirable  to  have  body-wave  phases  on  the  theoretical  seismograms; 
this  has  been  one  major  thrust  of  our  most  recent  work;  another 
has  been  the  improvement  in  efficiency,  control  of  accuracy  and 
increased  power  and  flexibility  of  the  computational  algorithms 
and  computer  programs.  We  have  reported  our  ability  to  generate 
theoretical  seismograms  with  body  and  surface  waves, on  the  same 
record  using  up  to  twenty-one  modes.  Multimode  theoretical 
seismograms  containing  body  wave  phases  have  been  generated  in 
the  past  (Sato,  Usami  and  Landisman,  1968) ; however,  these  time 
series  were  limited  to  ultralong  periods.  The  important  point 
of  our  recent  results  is  that,  owing  to  the  efficiency  of  our 
new  algorithms,  we  have  been  successful  in  extending  the  period 
content  of  the  theoretical  seismograms  through  the  range  covered 
by  the  WWSSNLP  instruments.  Thus,  we  can  generate  the  theoretical 
time  series  which,  for  the  first  time  we  believe,  permits  us  to 
compare  theoretical  seismograms  directly  with  the  entire  records 
obtained  at  the  WWSSN  and  the  HGLP  installations  down  to  a period 
of  ten  seconds.  This  requires  that  90-100  modes  be  used  in  the  multi- 
mode  synthesis.  The  necessary  frequency-domain  information  is 
obtained  in  a single,  relatively  inexpensive  computer  run;  time 
series,  for  any  source  specification,  are  then  obtained  with  a 
single  run  of  a second  program  (Liao,  Schwab  and  Mantovani  1977; 
preprint  appended) . 
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In  the  case  of  Rayleigh  waves,  sphericity  and  gravity  intro- 
duce complications  not  encountered  with  Love  waves.  For  Love 
waves,  optimized  flat-structure  programs  can  be  used  since  the 
sphericity  of  the  real  earth  can  be  handled  by  exact  transformation, 
and  since  gravity  does  not  affect  Love  waves.  At  present,  accurate 
Rayleigh-wave  computations  for  a real  earth  require  that  sphericity 
and  gravity  be  handled  directly  in  the  computational  algorithms 
and  programming;  hence,  our  first  task  here  was  to  improve  the 
accuracy  characteristics  of  direct  Rayleigh-,  or  spheroidal-wave 
computations.  The  difficulty,  expense,  and  time  required  for  this 
analysis  showed  us  why  this  information  was  not  in  the  literature, 
even  though  the  basic  algorithm  is  available  (Alterman,  Jarosch 
and  Pekeris,  1959).  The  results  of  our  work  (Schwab  et  al,  1977) 
are  appended.  They  include:  (1)  an  improved  and  simplified  com- 
putational algorithm  for  the  computation  of  the  phase  velocity  of  Rayleigh 
waves.  (2)  What  we  believe  to  be  the  first  direct  method  for 
computing  the  group  velocity,  i.e.,  without  appeal  to  variational 
methods;  (3)  numerical  analyses  and  examples  of  numerical  dif- 
ficulties encountered  in  this  type  of  computation;  (4)  detailed 
analysis  of  the  efficiency  of  our  optimized  algorithm  and  pro- 
gramming. Even  in  this  optimized  form,  these  direct  Rayleigh- 
wave  computations  for  a spherical,  gravitating  earth,  are  about 
six  times  slower  than  the  comparable  Love-wave  computations; 
whereas  we  know  that  computations  for  Rayleigh  waves  on  a flat, 
non-gravitating  structure  are  only  twice  as  slow.  Our  conclusion. 
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from  the  Rayleigh  wave  work  done  here,  i?  that  new  algorithms 
are  required  before  the  efficiency  of  the  Rayleigh-wave  com- 
putations (including  sphericity  and  gravity)  will  approach 

a level  justifying  computation  of  the  associated  multimode, 

i 

i 

} theoretical  seismograms  down  to  the  desired  period  of  ten 

I . 

seconds.  We  have  developed  a sufficiently  improved  algorithm  | 

(Schwab,  1977) , but  the  necessary  numerical  verification  and 
i t testing  are  only  in  their  initial  stages. 

I 

I References  to  our  own  work  on  theoretical  seismograms  for 

' Sa,  Lg  and  Sn  include:  Knopoff,  Schwab  and  Kausel  (1973); 

Knopoff  et  al  (1974) ; Nakanishi,  Schwab  and  Kausel  (1976) ; 

Kausel,  Schwab  and  Mantovani  (1977);  Mantovani  et  al  (1977). 

References  to  our  work  on  theoretical  seismograms  which  in- 
clude both  body  and  surface  waves  are:  Nakanishi,  Schwab 
and  Kausel  (1976) ; Kausel,  Schwab  and  Mantovani  (1977) ; 

! Mantovani  et  al  (1976);  Mantovani  (1977a, b);  Liao,  Schwab 

and  Mantovani  (1977) . 
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TECHNICAL  IMPORT 


I.  Purposes 

-f'  One  purpose  of  this  project  was  to  determine  the  regional  variations 
of  the  crust  and  upper  mantle  in  the  regions  under  investigation 
using  the  properties  of  surface  wave  dispersion  in  the  Arctic 
region  and  the  Eurasian  continental  area.  The  regions  of  high 
seismicity  within  and  around  the  region,]plus  the  dense  set  of 
WWSSN  stations  located  around  the  perimeter  of  the  region,  gave 
assurance  that  we  could  acquire  sufficient  single-station  phase 
velocity  data  for  the  application  of  regionalization  procedures. 

A second  purpose  was  to  develop  computational  techniques 
for  calculating  complete,  accurate  theoretical  seismograms  for 
both  earthquake  and  explosion  sources  within  the  area  under  in- 
vestigation, and  which  could  be  compared  with  those  recorded  at 
the  WWSSN  and  HGLP  stations  at  the  edge  of  this  region. 

II.  Review  of  scientific  background 

Regionalization.  The  single-station  surface  wave  method 
allows  all  stations  to  be  located  at  the  edge  of  the  region 
under  investigation;  this  technique  is  ideal  for  the  study  of 
Eurasia.  Knopoff  and  Schwab  (1968)  extended  the  description  of 
the  single  station  method  (Brune,  Nafe  and  Oliver,  1960)  to  take  into 
account  the  frequency  dependence  of  the  apparent  initial  phase 
of  the  source.  The  success  of  the  single-station  method 
depends  on  the  calculation  of  the  initial  phase.  For  dislocation 
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sources  in  layered  media,  the  initial  phase  can  be  determined 
by  a method  due  to  Ben  Menahem  and  Harkrider  (1964)  and 
Harkrider  (1964,  1970)  who  obtained  the  surface  wave  response  to 
double  couple  sources.  The  far-field  response  to  displace- 
ment-dislocation faulting  is  equivalent  to  that  from  a point- 
source,  double-couple  in  an  unfaulted  medium  (Burridge  and 
Knopoff , 1964) . Thus  if  the  source  focal  mechanism  is  known 
from  fault-plane,  surface  wave  amplitudes,  or  other  methods, 
the  initial  phase  can  be  determined  for  small  or  moderate 
sized  earthquakes.  By  means  of  transformation  techniques 
(Biswas  and  Knopoff,  1970)  or  empirical  correction  methods 
(Bolt  and  Dorman,  1961) , it  is  possible  to  provide  corrections 
for  sphericity. 

The  first  inversions  of  surface  wave  dispersion  data  to 
give  upper  mantle  structure  were  carried  out  by  Dorman  and 
Ewing  (1962),  Brune  and  Dorman  (1964)  and  Knopoff,  Mueller  and 
Pilant  (1966) . These  papers  were  concerned  with  obtaining  a 
single  structure  which  fit  the  experimental  data.  Knopoff 
(1961,  1962)  pointed  out  that  the  inversion  of  noise-free 
data  is  not  unique.  Subsequent  efforts  were  mainly  concerned 
with  finding  the  set  of  structural  models  which  fit  the  data  to 
within  the  experimental  accuracy  (Keilis-Borok  and  Yanovskaya, 
1967;  Press,  1968,  1969).  The  inversion  of  noisy  data  enlarges 
the  span  of  non-unique,  acceptable  models  over  that  for  noise- 
free  data. 

Two  methods  of  inversion  are  available.  In  the  first  of  the 


A 
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techniques,  perturbation  methods  are  applied  in  a linearized 
version  of  the  problem.  The  problems  of  inconsistent  or 
correlated  data  are  accounted  for  by  smoothing  in  the  data- 
space,  usually  by  a least-squares  procedure.  The  problems  of  over- 
parameterization of  the  model  space  are  taken  care  of  by  smoothing 
in  the  model  space,  by  a procedure  usually  described  as  "deltaness". 
Further  reductions  of  the  number  of  model  parameters  can  stabilize 
the  inversions.  Examples  of  linearized  inversions  are 

given  by  Backus  and  Gilbert  (1968,  1970),  Knopoff  & Jackson  (1972) 
Jackson  (1972) , and  Wiggins  (1972) . 

The  full  non-linear  inverse  problem  in  a multidimensional 
parameter  space  of  high  order  has  been  attacked  by  Monte  Carlo 
methods  (Press,  1968,  1969).  The  inefficiency  of  Monte  Carlo 
methods  can  be  minimized  at  the  sacrifice  of  reducing  the 
dimensionality  of  the  parameterization  to  a value  less  than  about 
seven  by  exploring  a neighborhood  of  acceptable  solutions  for  other 
acceptable  solutions.  Such  a program  has  been  called  Hedgehog 
and  has  been  much  used  in  this  laboratory  (e.g.  Biswas  and 
Knopoff,  1974  ; Knopoff  and  Schlue,  1972  ; etc.);  copies  of 
this  program  exist  in  Moscow,  Bologna,  Edmonton,  Bari,  Cambridge  (UK),  Paris, 
etc.  Both  linearized  and  Hedgehog  methods  have  been  used  in  our 
inversions. 
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The  single-station,  surface  wave  regionalization  of  an 
area  is  based  on  the  assxomption  that  the  surface  waves  travel 
paths,  and  the  assumption  that  the  phase  travel  time  from  epi- 
center to  station  is  the  sum  of  the  travel  times  through  the 
homogeneous  subdivisions  of  the  laterally  heterogeneous  region 
(Knopoff , 1969) . This  permits  the  travel-time  at  each  frequency 
to  be  expressed  as  a system  of  inhomogeneous  equations; 
for  each  path,  the  total  travel  time  is  the  inhomogeneous  term, 
the  distances  through  the  subdivisions  are  the  coefficients,  and 
the  slownesses  in  the  subdivisions  are  the  unknowns.  We  have 
pointed  out  (Leeds  et  al . 1974)  that  the  solution  of  this  system 
of  equations  can  yield  the  experimental  phase  slownesses 
associated  with  each  of  the  subdivisions  only  if  we  make  the 
assumption  that  the  errors  in  each  region  are  uncorrelated. 

Since  these  errors  are  not  uncorrelated,  we  must  consider  the 
model  parameters  as  the  primary  unknowns  in  the  inversion  and 
derive  the  phase  slownesses  in  the  pure  regions  therefrom. 

The  first  application  of  the  use  of  single-station  phase 
delay  measurements  in  the  pursuit  of  regionalization  of  a large 
inhomogeneous  area  was  performed  by  Kausel,  Leeds  and  Knopoff 
(1974) , Leeds,  Kausel  and  Knopoff  (1974)  and  Leeds  (1975) . In 
the  above  work,  phase  delays  along  long  paths  across  the  Pacific 
area  were  found  to  vary  systematically  with  distance  from  the 
East  Pacific  Rise.  This  persuaded  us  that  the  appropriate 
geographic  regionalization  was  plausibly  based  on  a basement 
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i.e.  magnetic,  age  of  sea-floor  spreading.  Our  geographic 
provinces  then  took  the  form  of  broad  magnetic  age  stripes. 

An  inversion  using  this  regionalization  showed  that  although  the 
number  of  criss-crossing  paths  was  plentiful,  the  number  of 
degrees  of  freedom  in  the  data  was  remarkably  small.  Furthermore, 
we  found  that  the  data  set  did  not  permit  one  to  obtain  detailed 
information  concerning  the  bottom  of  the  low-velocity  channel.  We  did 
find  that  the  lithosphere  increased  in  thickness  inonotonically  with  age; a 
the  ridge  crest,  the  lid  of  the  channel  has  almost  zero  thickness, 
while  in  the  oldest  ocean  this  lid  is  about  100  km  thick.  The  lid  thick- 
nesses are  consistent  with  a geochemical  model  in  which  the  lid- 
channel  interface  is  at  the  solidus  for  wet  peridotite. 

Because  they  have  not  developed  good  long-period  instruments, 

Soviet  seismologists  have  not  been  able  to  focus  attention  on 
upper  mantle  studies.  Also,  their  extensive  program  in  deep  seismic 
sounding  has  focused  interest  upon  crustal  studies.  The  Soviet 
surface  wave  work  has  been  limited  to  short-period  investigations  — 
usually  less  than  40  seconds  --  and  has  been  concerned  mainly  with 
determining  crustal  properties.  References  to  Soviet  surface-wave 
work  include  Arkhangelskaya  (1960) , Savarensky  and  Ragimov  (1958, 

1959),  Savarensky,  Solov'eva  and  Shechkov  (1959)  and  Savarensky 
and  Sikharulidze  (1959),  Popov  (1960),  Shechkov  (1961,  1964,  1970), 
Savarensky  and  Shechkov  (1961) , Shechkov  and  Solov'eva  (1961) , 

Savarensky  and  Peshkov,  1968;  Sikharulidze  and  Makharadze,  1968; 
Savarensky  et  al . 1969. 


since  appropriate  Soviet-based  seismic  data  are  not  avail- 
able to  us,  our  own  work  has  thus  far  focused  on  the  estimation 
of  upper  mantle  properties  in  Eurasia  through  the  use  of  the 
only  reasonable  tools  available  to  us,  namely  phase  delays  of 
surface  waves  on  long  paths,  criss-crossing  the  region  of 
interest.  This, plus  a plausible  model  of  regionalization  — 
based  in  part  on  the  observations  in  Knopoff  (1972)  for 
continents  — , permits  an  attack  on  the  problem  of  regionalization 
of  the  upper  mantle  of  Eurasia. 

Theoretical  seismograms.  Relative  to  the  discrimination 
problem,  probably  the  most  important  feature  in  the  calculation 
of  theoretical  seismograms  which  requires  improvement  over 
previously  existing  systems  for  such  computations  is  the  capa- 
bility of  extending  both  body-  and  surface-wave  portions  of  the 
theoretical  computer  seismograms  to  short  periods.  In  this 
context,  by  "short-period"  we  refer  to  the  period  range  covered 
by  the  long-period  instruments  of  the  WWSSN  installations.  In 
our  formulation,  the  successful  accomplishment  of  this  task  is 
dependent  upon  improved  techniques  for  obtaining  multimode 
dispersion-attenuation  information  for  reasonably  realistic 
models  of  the  earth,  i.e.  spherical,  radially  heterogeneous, 
anelastic  models. 

Optimization  of  this  type  has  been  one  of  the  main  interests 
in  our  laboratory  for  several  years.  The  results  of  our  early 
work,  based  on  the  Thomson  (1950) -Haskell  (1953)  technique  and 
on  Knopoff 's  (1964)  method  for  treating  flat-layered  structures. 


are  reported  by  Schwab  (1970)  and  Schwab  and  Knopoff 
(1970;  1971;  1972;  1973).  The  results  of  our  work  on 
spherical-to-f lat-structure  transformation  techniques,  which 
permit  the  use  of  flat-structure  programs  in  dispersion- 
attenuation  computations  with  Love  waves  for  spherical  models 
of  the  earth,  are  given  by  Biswas  and  Knopoff  (1970) , Schwab 
and  Knopoff  (1971;  1972;  1973) , and  Kausel  and  Schwab  (1973) . 

In  this  last  reference,  we  have  also  given  an  outline  of  the 
approach  we  have  adopted  to  handle  the  synthesis  of  multimode 
seismograms  once  the  dispersion,  attenuation,  source,  and 
excitation  functions  have  been  specified.  The  entire 
theoretical  seismogram  for  a dislocation  source  in  a spherical 
earth  can  be  expressed  as  a simple  sum  of  normal  mode  contri- 
butions (Saito,1967;  Takeuchi  and  Saito,  1972).  We  first  applied 
our  scheme  for  generating  theoretical  seismograms  to  the  inter- 
pretation of  the  seismic  phase  Lg  (Knopoff,  Schwab  and  Kausel, 

1973;  Knopoff  et  al.  1974)  This  phase  is  a multimode  interference 
pheonomenon  which  belongs  to  the  surface  wave  portion  of  the  seismo- 
gram. In  this  report  we  indicate  that  we  have  developed  a system 
for  synthesis  where  both  body  and  surface  waves  appear  on  the 
same  seismogram  for  realistic  models  for  the  earth,  and  where  the 
period  range  spans  that  covered  by  the  long-period 
instruments  of  the  WWSSN.  Earlier  work  of  this  type,  which  was 
performed  with  simplified  models  of  the  earth,  is  summarized  by 
Alterman  and  Loewenthal  (1972) . Sat6,  Usami  and  Landisman  (1968) 
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describe  the  computation  of  complete  theoretical  seismograms 
for  realistic  models  of  the  earth.  However,  their  results  are 
limited  to  ultralong  periods. 

III.  Objectives,  methods  and  results. 

Data  collection.  Accurate  surface  wave  phase  delays 
across  the  area  with  single-station  method  require  that  we 
know  the  focal  mechanism  and  the  depth  of  focus  in  order  to 
obtain  the  corrections  due  to  the  apparent  initial  phase.  We 
have  confined  our  attention  to  the  measurement  of  the  phase 
velocity  dispersion  of  the  fundamental  mode.  Frez  and  Schwab 
(1976)  have  computed  the  importance  of  the  structural  parameters 
on  the  determination  of  the  initial  phase. 

Long  period  records  from  the  47  WWSSN  stations  which 
border  the  region  of  interest  have  been  used  in  the  study. 

The  locations  of  these  stations  are  shown  in  Figure  1.  We  have 
also  obtained  several  seismograms  (through  World  Data  Center  B) 
of  Soviet  records  made  in  Central  Asia  on  experimental  long- 
period  instruments.  However,  the  precision  of  these  latter 
recordings  plus  an  uncertain  calibration  impulse  response  has 
not  permitted  us  to  use  them. 

An  example  of  the  intermediate  magnitude  seismicity  of 
the  area  is  given  in  Figure  2.  Epicenters  are  plotted  for 
earthquakes  which  occurred  during  the  interval  from 
February,  1963  to  February,  1967  and  having  magnitudes  between 
5.9  and  6.6.  In  addition  to  the  epicenters  shown  in  Figure  2, 
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FIG.  1.  Locations  of  WWSSN  and  HGLP  stations. 

FIG.  2.  Seismicity  of  the  Eurasian  region.  The 

regions  of  high  seismicity  along  the  eastern 
border  of  the  Kamchatka  peninsula  and  along 
the  Aleutian  arc  also  provide  useful  events 
for  the  study. 
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ISEFUL  EPICENTERS 
i.9<  MAGNITUDE  <6.6 
EB.  1963  - FEB.  1967 
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there  are  regions  of  high  seismicity  along  the  eastern  border 
of  the  Kamchatka  peninsula  and  along  the  Aleutian  arc.  The 
choice  of  this  range  of  magnitudes  is  governed  by  two  considerations. 
First,  experience  has  shown  that  good  long-period  surface  wave 
information  requires  events  with  magnitudes  above  a certain  value; 
of  course,  a shock  which  is  so  large  as  to  send  the  instrument 
off-scale  is  useless  for  our  purposes.  Second,  the  application 
of  the  single-stations  method  requires  knowledge  of  the  focal 
mechanism.  We  must  therefore  use  events  large  enough  to  allow 
us  to  obtain  an  accurate  fault  plane  solution  for  each  event  we 
select  for  processing. 

Since  the  set  of  stations  around  the  area  to  be  studied 
is  dense  as  are  the  regions  of  high,  intermediate-magnitude 
seismicity  located  within  and  around  the  area,  there  has  been 
no  problem  in  obtaining  sufficient  data  for  the  project.  It  is 
interesting  to  note  that  the  area  is  almost  completely  encircled 
by  either  stations  or  epicenters  or  both.  The  limits  of  the 
area,  which  we  have  covered  with  a sufficiently  dense  set  of 
paths  from  earthquakes  to  epicenters,  are  shown  in  Figure  3 
by  the  solid  lines.  The  shaded  regions  are  those  of  high,  large- 
magnitude  seismicity. 

Recordings  from  fifteen  events  were  processed.  The  list 
of  events  is  given  in  Table  I,  with  their  USCGS-NOAA  specifi- 
cations. The  focal  parameters  and  other  data  are  listed  in 
Table  II.  We  have  constructed  fault  plane  solutions  for  11 


Table  I 


Events  used  in  phase  velocity  determinations 
(USCGS  specifications) 


Date 

Time 

Lat. 

Long. 

h (km) 

M 

1. 

Red  Sea 

Mar. 

31, 

1969 

07:15:54.4 

27.7°N 

34.0°E 

33 

6.0 

2. 

Hsingtai 

Mar. 

7, 

1966 

21:29:17.4 

37.3°N 

114. 9°E 

33 

6.0 

3. 

Kamchatka 

Dec. 

26, 

1964 

14:30:29.1 

51. 8°N 

156. 8°E 

136 

5.7 

4. 

Lena  River 

Aug. 

25, 

1964 

13:47:20.6 

78.2°N 

126. 6°E 

50 

6.1 

5. 

Tashkent 

Apr. 

25, 

1966 

23:22:49.3 

41.3°N 

69.2°E 

8 

5. 

6. 

Lop  Nor 

Oct. 

14, 

1970 

07:29:58.6 

40.9°N 

89.4°E 

0 

— 

7. 

Yunnan 

Feb. 

13, 

1966 

10:44:41.3 

26.1°K 

103. 2°E 

33 

5.7 

8. 

Hindu  Kush  I 

June 

6, 

1966 

07:46:16.1 

36.4°N 

71.1°E 

221 

6.2 

9. 

Hindu  Kush  II 

Dec. 

28, 

1974 

12:11:43.7 

35.1°N 

72.9°E 

22 

6.0 

10. 

East.  Aleut. 

Feb. 

6, 

1965 

01:40:33.2 

53.2°N 

161. 9°W 

33 

6.4 

6 

11. 

East.  Aleut. 

Feb. 

6, 

1965 

16:50:23.6 

53. 3°N 

161. 8°W 

33 

6.1 

6 

12. 

West.  Aleut. 

Feb. 

5, 

1965 

09:32:09.3 

52. 3°N 

174. 3°E 

41 

5.9 

6 

13. 

West.  Aleut. 

Feb. 

6, 

1965 

04:02:53 

52.1°N 

175. 7°E 

35 

5.9 
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Table  II 

Focal  Parameters  for  Events  Studied 


Fault  Plane 


ent  No. 

<{> 

6 

X 

h (km) 

Solutions 

Paths 

1. 

133° 

66° 

242° 

15 

Fig.  4 

Fig.  14 

2. 

122° 

82° 

0° 

14.5 

5 

15 

3. 

210° 

90° 

90° 

136. 

6 

16 

4. 

165° 

58°W 

260°^'^ 

11.5^'^ 

7 

17 

5. 

305°^ 

70°^ 

90°^ 

8. 

22 

6. 

explosion 

0. 

22 

7. 

190° 

70° 

279° 

33. 

8 

18 

8. 

77° 

50° 

76° 

221. 

9 

18 

9. 

122° 

4 0° 

61° 

22. 

10 

18 

10. -11. 

(twin 

earthquakes) 

11a, b 

19 

12.-13. 

(twin 

earthquakes) 

12 

20 

Depth 

obtained 

from  Rayleigh  wave 

spectra . 

Rayleigh  wave 

spectra 

for  event  2 

shown  in 

Fig.  13 

Sykes 

(1967)  gives  fault  planes  for  this  event  as 

(j)  = 

4°, 

6 = 58°W 

<J)  = 

338°, 

6 = 54°E 

19 
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I ‘ 

i 

i 

of  these  cases  (including  two  pairs  of  twin  earthquakes) . The 
\ list  of  figures  of  these  fault  plane  solutions  is  given  in 

Table  II.  In  one  case  (No.  5)  we  have  used  the  fault  plane 
solution  of  Zakharova  et  al  (1971) . In  one  case  (No.  6)  we 

I ' have  assumed  the  source  that  of  an  explosion.  In  two  cases 

I 

(Nos.  1,  2)  we  have  used  Rayleigh  wave  spectra  to  help  deter- 

[ mine  the  focal  parameters;  in  these  cases,  we  were  able  to 

r 

j obtain  the  depth  of  focus  more  accurately  than  by  conventional 

^ methods.  In  the  case  of  No.  2 the  spectra  (Fig.  13)  also  helped 

refine  the  fault  plane  parameters  and  permit  us  to  modify  Sykes' 
(1967)  values.  The  methods  used  to  determine  fault-plane 
parameters  from  surface  wave  spectra  are  given  by  Ben  Menahem 
and  Toksoz  (1963) . 
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FIG.  4.  Fault  plane  solution  for  event  4 (Red  Sea) 

occurring  at  07:15:54.4  GMT,  March  31,  1969. 

FIG.  5.  Fault  plane  solution  for  event  1 (Hsingtai) 
occuring  at  21:29:17.4  GMT,  March  7,  1966. 

FIG.  6.  Fault  plane  solution  for  event  3 (Kamchatka) 

occuring  at  14:30:29.1  GMT,  December  26,  1964. 

FIG.  7.  Fault-plane  solution  for  event  2 (Lena  River) 
occurring  at  13:47:20.6  GMT,  August  25,  1964. 

The  solution  given  by  Sykes  (1967)  is  indicated 
by  dotted  lines. 

FIG.  8.  Fault  plane  solution  for  event  7 (Yunnan) 

occurring  at  10:44:41.3  GMT,  February  13,  1966. 

FIG.  9.  Fault  plane  solution  for  event  8 (Hindu  Kush) 
occurring  at  07:46:16.1  GMT,  June  6,  1966. 

FIG.  10.  Fault  plane  solution  for  event  9 (Hindu  Kush) 
occurring  at  12:11:43.7  GMT,  December  28,  1974. 

FIG.  11a.  Fault  plane  solutions  for  events  10  and  11 
(Eastern  Aleutians)  occurring  at  01:40:33.2 
(Fig.  11a)  and  16:50:29  GMT  (Fig.  11b),  February  6, 


1965. 
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FIG.  12.  Fault  plane  solution  for  events  12  and  13 

(Western  Aleutians)  occurring  at  09:32:09.3  GMT, 
February  5,  1965  and  04:02:53  GMT,  February  6, 
1965. 


FIG.  13.  Rayleigh-wave  amplitude  distributions  for 

event  2 (Lena  River)  occurring  at  13:47:20.6 
GMT,  August  25,  1964.  The  central  set  of 
radiation  patterns  are  the  results  of  theor- 
etical computations  based  on  the  fault  plane 
solution  given  in  the  text.  The  other  four 
radiation  patterns  depict  the  experimental 
results . 


A 


19H 
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Two  points  should  be  noted  concerning  the  accuracy  of  the 
digitizations  of  the  recorded  events.  First,  most  of  the  event 
records  we  have  used  are  about  as  large  as  they  could  be  without 
going  off  scale.  This  has  necessitated  a change  in  our  data 
processing  techniques  which  should  be  noted  for  the  information 
of  others  involved  in  this  type  of  work. 

In  the  past,  our  standard  procedure,  when  working  with 
smaller-amplitude  recordings,  has  been  to  digitize  from  copies 
of  35  mm  microfilm  records  of  the  WWSSN  seismograms  made  with  a 
standard  microfilm  reader-printer  (Itek  18.24  Reader-Printer). 
Tests  which  compare  the  phase  velocity  results  obtained  from 
full-size  record  copies  provided  by  NOAA  with  the  results 
obtained  from  our  microfilm  copies  show  that  distortion  in  the 
copying  process  is  of  concern  when  working  with  large-amplitude 
recordings  such  as  those  employed  in  the  present  study.  We  en- 
courage the  use  of  full-size  record  copies  of  large  events 
obtained  directly  from  NOAA. 

The  second  point  concerns  the  fact  that  the  direction  of 
swing  of  the  galvanometer  may  not  be  parallel  to  the  axis  of 
the  recording  drum.  Although  James  and  Linde  (1971)  term  this 
phenomenon  "a  source  of  major  error  in  digital  analysis  of  WWSSN 
seismograms",  our  tests  show  the  effect  to  be  negligible  on  phase 
travel  times  computed  using  the  single-station  method  for  epicenter- 
station  separations  of  a few  thousand  kilometers.  In  the  case  of 
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the  poorest  galvanometer  alignment  we  encountered,  about  triple 

I the  normal  ramp  slope  of  0.3°,  we  found  only  negligible  differences 

I ' between  the  phase  velocity  curve  obtained  with  the  correct 

digitization  base  line,  and  the  curve  obtained  with  the  normal 
ramp  as  a base  line. 

The  paths  over  which  the  phase  delays  have  been  measured  from 
each  of  the  seismic  events  are  shown  in  Figs.  14-20  and  are  sum- 
marized in  Table  II.  Sample  phase  velocity  results  are  given  in 

I Figure  21,  which  illustrates  the  variation  in  dispersion  for  dif- 

1 

ferent  propagation  paths.  We  have  obtained  phase  velocity  data  for 
most  paths  over  a period  range  extending  from  about  30  or  38  sec. 
in  most  cases,  to  as  long  as  357  sec  in  a few  rare  cases.  The 
instrumental  response  at  these  longest  periods  is  unreliable;  if 
the  measured  values  of  phase  velocity  and  its  first  derivative 
were  in  significant  disagreement  with  values  from  the  free  mode 
spectrum,  then  these  long  period  values  were  rejected.  For  our 
present  inversions,  we  have  only  used  periods  up  to  250  sec.  The 
specific  period  ranges,  which  we  have  used  in  the  inversions  are 
shown  in  Table  3. 

The  phase  velocities  from  five  earthquakes  and  one  nuclear 
explosion  for  the  32  paths  crossing  Eurasia  (Figure  22)  sort 
themselves  into  two  groups  (Figures  23  and  24) . The  paths  with 

j higher  phase  velocities  are  generally  those  that  cross  the  stable 

i platforms  and  shields  (such  as  paths  from  the  Hsingtai  earthquake 

I 

to  Scandinavian  (KEV)  and  German  (STU)  stations;  typical  of  the 

i 

I lower-velocity  group  is  the  phase  velocity  on  the  paths  from  the 

Red  Sea  earthquake  to  the  southern  Asiastic  stations  (SHL,  MAN,  etc.)). 
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FIG. 


14.  Epicenter-to-station  paths  processed  for 
event  4 (Red  Sea)  occurring  at  07:15:54.4  GMT, 

March  31,  1969. 

15.  Epicenter-to-station  paths  processed  for 
event  1 (Hsingtai)  occurring  at  21:29:74.4  GMT, 

March  7,  1966. 

16.  Epicenter-to-station  paths  processed  for 

event  3 (Kamchatka)  occurring  at  14:30:29.1  GMT, 

December  26,  1964.  j 


i 

17.  Epicenter-  to-station  paths  processed  for  i 

event  2 (Lena  River)  occurring  at  13:47:20.6  GMT,  | 

August  25,  1964.  ! 

18.  Epicenter-to-station  paths  processed  for  event 
7 (Yunnan)  occurring  at  10:44:41.3  GMT, 

February  13,  1966,  and  events  8 and  9 

(Hindu  Kush)  occurring  at  07:46:16.1  GMT, 

June  6,  1966  and  12:11:43.7  GMT,  December  28, 


1974. 
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FIG. 


FIG. 


19.  Epicenter-to-station  paths  processed  for 
events  10  and  11  (Eastern  Aleutians)  occurring 
at  01:40:33.2  (dashed  lines)  and  16:50.29  GMT 
(solid  lines),  February  6,  1965. 

20.  Epicenter-to-station  paths  processed  for 
events  12  and  13  (Western  Aleutians)  occurring 
at  09:32:09.3  GMT,  February  5,  1965  (dashed 
line)  and  at  04:02:53  GMT,  February  6,  1965 
(solid  line) . 

21.  Sample  phase  velocity  results  for  the  paths 
from  the  Hsingtai  earthquake  (1966)  to  STU 
and  from  the  Lena  River  earthquake  (1964)  to 
HOW.  The  paths  are  shown  in  Figures  15  and 
16. 
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PERIOD  (seconds) 


TABLE  3 

PERIOD  RANGE  (sec) 

PATH 

30 

38 

50  69  100  119 

RED  SEA-ANP 

X — 

HSINGTAI-ATU 

-KEV 

-LAH 

-MSH 

-NDI 

-POO 

-STU 

-TAB 


KAMCHATKA- CHG 


LENA  RIVR-ATU 


TASHKENT  -NDI 


LOP  NOR  -KBL 
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YUNAN 


HINDU  KUSH  I 


HINDU  KUSH  II 


TABLE  3 

PERIOD  RANGE  (sec) 

PATH  30  38  50  69  100  119  139  167  192  208  227  250 

ANP  X X 

BAG  X ^X 

MSH  X ^X 

SEO  X X 


•SHK  X X 

ANP  X X 

HKC  X X 

ANP  X X 
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FIG.  22.  Propagation  paths  across  the  Eurasian  continent 
from  five  earthquakes  and  one  nuclear  explosion 

FIG.  23.  All  Eurasian  phase  velocities  can  be  sorted  into 
two  groups  (shaded  areas) , except  for  phase 
velocities  POO-2  and  SHL-1,  which  fall  between 
these  two  groups.  Phase  velocities  for 
"standard"  shield  (FLO-GOL) , younger  stable 
regions  (SHA-LUB)  and  rift  zones  (TUC-BOZ) 
are  shown  for  comparison  (Biswas  and  Knopoff, 
1974)  . The  global  average  phase  velocities 
obtained  from  free-mode  observations  are  also 
shown  (F.M.)  (Gilbert  and  Dziewonski,  1975). 

FIG.  24.  Propagation  paths  corresponding  to  the  two 

phase  velocity  groups  in  Fig.  21.  The  solid 
line  indicates  a path  with  higher  phase  velocity; 
the  dashed  line  indicates  a path  with  lower 
phase  velocity. 
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Phase  velocities  for  typical  shield  regions  and  young  stable  con- 
tinental regions  (Biswas  and  Knopoff,  1974)  are  shown  for  com- 
parison. The  incompatibility  of  most  of  the  phase  velocity 
observations  for  Eurasian  paths  with  the  "standard"  curves  for 
homogeneous  regions  testifies  to  the  inhomogeneous  nature  of  the 
Eurasian  region. 

Regionalization  i 

1 

The  first  step  in  the  structural  analysis  requires  the  j 

division  of  the  area  into  subregions;  each  of  these  is  assumed 

I 

to  be  laterally  homogeneous.  In  our  first  attempt  we  have  sub-  ' 

divided  the  Eurasian  area  into  six  broad  regions  (Fig.  25) . These  are: 

j 

1.  Ancient  PreCambrian  Shields 

2.  Stable  Platforms 

3.  The  Himalayan-Alpide  Mountain  Belt 

4.  The  Tibetan  Plateau 

5.  The  Sinkiang-Mongolian  Seismic  Zone 

6.  The  Chinese  "Stable"  Region  i 

This  choice  of  regionalization  is  largely  governed  by  an  attempt 

to  define  a small  number  of  discrete  geographically  homogeneous 
provinces.  In  this  choice  we  have  been  guided  by  tectonic  maps 
of  Eurasia  (Khain  and  Muratov,  1969)  as  well  as  by  seismicity  and 
sparse  heat  flow  information (Lubimova  and  Polyak,  1969).  Our  first 
impulse  was  to  regionalize  the  stable  parts  of  Eurasia  according  to  j 

the  bimodal  division  of  these  regions  into  stable  ancient  shields  j 

and  stable  younger  continents  (Knopoff,  1972).  However,  I 

information  about  basement  ages  | 

i 

i 

I 

J 
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of  the  Eurasian  sediment-covered  platforms  is  difficult  to  obtain, 
so  we  have  taken  the  simple  expediment  of  dividing  the  stable 
parts  of  the  continent  into  the  geologically  identif icable  shields 
and  a second  region  which  includes  the  rest  of  the  seismically  stable  land 
mass;  these  form  our  regions  1 and  2.  In  selecting  region  2 as  a 
homogeneous  region  we  impose  no  a priori  condition  that  it  be  either 
of  the  north-central  U.S.  or  the  Gulf  Coast  (U.S.)  types  of 
sediment-covered  stable  regions:  the  first  of  these  regions  in  the 
U.S.  has  been  found  to  have  an  upper  mantle  with  shield  character- 
istics while  the  second  has  an  upper  mantle  with  a strong  low  velocity 
channel,  and  is  the  prototype  younger  stable  region. 

The  next  three  regions  are  characterized  by  their  high 
tectonic  activity.  We  have  chosen  to  identify  the  mountainous 
collision  zone  between  the  Asian  and  Indian  plates  as  a single 
province;  this  is  undoubtedly  incorrect  in  detail,  but  is  valid 
enough  in  view  of  our  inability  to  provide  detailed  resolution  of 
smaller  regions  by  our  use  of  long  paths. 

It  should  be  noted  that  an  increase  in  the  density  of  path 
samples  does  not  necessarily  improve  the  resolution  of  structure 
of  a small  geographic  area.  Consider  the  limiting  case  of  a 
region  which  is  vanishingly  small  in  area.  An  increase  in  the 
number  of  paths  which  transect  this  region  does  not  improve  our 
ability  to  estimate  the  cross-section  under  it,  because  the 
travel-time  spent  by  the  ray  in  this  region  is  a vanishingly 
small  part  of  the  total  delay.  An  increase  in  the  number  of  paths 
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in  a small  area  only  increases  the  redundancy  of  the  data.  For  this 
reason  we  had  some  hesitation  regarding  the  possibility  of  determin- 
ing the  structure  under  a region  as  "small"  as  the  Tibetan  plateau. 
The  Tibetan  plateau  is  indeed  small  in  comparison  with  the  dimen- 
sions of  the  Eurasian  continent,  and  herein  lies  one  of  the  dif- 
ficulties with  the  method  of  very  long  single-station  transverses 
across  the  entire  span  of  the  continent.  But  the  Tibetan  plateau 
contains  one  of  the  important  mysteries  of  modern  Plate  Tectonics, 
namely  the  reason  for  its  great  elevation  and  the  nature  of  the 
architectural  underpinning  that  holds  it  up,  so  we  resolved  to  try 
to  determine  its  structu  e, fully  anticipating  that  the  error  bias 
in  the  determinations  of  structural  parameters  might  be  large;  it 
became  our  region  4.  To  improve  resolution  here,  we  made  use  of 
shorter  paths  from  two  Hindu  Kush  and  one  Tashkent  earthquakes  plus 
a Lop  Nor  explosion  all  recorded  at  nearby  stations  to  increase  the 
fraction  of  the  travel- times  spent  in  the  Tibetan  region.  It  can  be 
seen  from  the  table  of  period  ranges  of  the  observations  (Table  3) 
that  the  inability  to  obtain  long-period  information  from  the 
observations  of  the  Tashkent  earthquake  and  the  Lop  Nor  nuclear  ex- 
plosion also  limits  our  ability  to  resolve  deeper  structure  under 
region  4.  The  complete  sampling  of  all  regions  by  the  path-phase 
delays  from  the  set  of  earthquakes  and  explosions  is  shown  in  Fig.  26, 
The  Sinkiang-Mongolia  seismic  zone  is  easily  identified  as  one 
of  the  provinces  in  our  regionalization.  We  have  chosen  to  identify 
the  relatively  stable  aseismic  block  of  Southeast  China 
as  an  additional  province.  The  Hsingtai  earthquakes  of  1968 


occurred  on  the  boundary  between  these  two  regions.  We  have 
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identified  Southeast  China  as  a separate  region  without  any  a priori 
assumption  that  it  has  either  a cross-section  similar  to  those  for 
the  ancient  shields  or  the  younger  stable  continents.  We  prefer  to 
allow  the  inversion  to  yield  this  determination.  Should  the  in- 
version indicate  that  this  region,  or  region  2 for  that  matter,  is 
similar  in  cross-section  to  that  for  some  other  region,  we  may  make 
the  assumption  that  these  regions  are  the  same  and  use  this  information 
to  perform  a further  simplified  inversion.  Although,  within  them- 
selves, these  regions  encompass  widely  differing  geologic  structures 
and  widely  varying  seismic  activity,  we  have  assumed  that  each  of 
these  regions  is  homogeneous  in  order  to  limit  the  number  of 
parameters  in  the  inversion. 

Three  small  regions  are  treated  specially.  For  the  South 
China  Sea,  we  have  assumed  the  structure  to  be  known,  and  to  be 
that  for  typical  marginal  seas.  The  phase  delays  for  this  region 
are  taken  from  two  single-station  phase  velocity  observations  in 
marginal  seas  obtained  by  Leeds  (1973)  and  from  two  phase  velocity 
determinations  by  the  two-station  method  across  the  Philippine  Sea. 
Phase  velocity  corrections  for  the  Sea  of  Okhotsk  were  obtained 
theoretically  from  a crustal  structure  given  by  Kosminskaya  et  al 
(1969)  and  derived  from  explosion  work,  in  which  the  Okhotsk 
depression  has  a 25  km  crustal  thickness;  we  have  used  an  oceanic 
mantle  below  this  crust;  These  phase  corrections  have  been  applied 
to  travel  times  for  those  paths  that  traverse  these  two  regions. 


26 


Finally,  a small  fraction  of  the  paths  cross  the  Yakutsk-Magadan 
region  of  northeastern  Siberia.  The  maps  show  this  mountainous 
region  to  be  tectonically  different  in  a significant  way  from 
the  shield  immediately  to  the  west  of  it.  Since  our  samples  of 
this  region  are  all  from  the  Kamchatka  earthquake,  and  are  all 
small  in  length,  we  have  decided,  without  justification,  to 
couple  this  region  to  our  other  mountainous  province,  region  3. 

We  do  not  expect  significant  differences  from  a region  3 structure 
for  this  part  of  Eastern  Eurasia  to  produce  major  changes  in  the 
inversions  since  the  total  fraction  of  path  length  in  this  region 
is  small. 

The  total  path  length  in  each  region  summed  over  all  event- 
paths  in  given  in  Table  4.  Since  the  phase  shifts  have  not  been 
obtained  over  a uniform  band  of  periods  for  all  paths,  in  the  last 
column  of  Table  4 we  have  also  reported  an  estimate  of  the 
number  of  samples  in  each  region  by  giving  the  product  of  sample 
path  length  by  number  of  period  estimates  of  phase  shift.  By 
either  method  of  estimation,  the  very  low  fraction  of  sampling 
in  regions  4 and  6 lead  us  to  expect  that  large  uncertainties  in 
the  structure  will  be  obtained  from  the  inversion  for  these 
regions . 

We  have  inverted  the  data  under  the  assumption  that  a simple 
ray  theory  for  surface  waves  applies,  that  is,  the  phase  shift  for 
a surface  wave  passing  through  a given  region  is  computed  as  though 
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TABLE  4 


Reqion 

Path 

Length (km) 

Percentage  of  Total 
Path  Length 

Weighted  Percentage 
Total  Path  Length 

1 

20970 

12.1 

13.0 

2 

53150 

30.6 

30.4 

3 

38411 

22.1 

22.5 

4 

15292 

8.8 

6.9 

5 

37994 

20.7 

20.7 

6 

9979 

5.7 

6.5 

TOTAL  173796  100.00  100.0 
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the  region  were  laterally  infinite  in  extent  and  uninfluenced  hy 


the  presence  of  neighboring  regions,  no  matter  how  close  the  great 
circle  path  approaches  the  regional  boundaries.  The  assumption 
that  diffraction  effects  are  unimportant  is  evidently  untenable 
but  provides  a basis  for  starting  an  inversion;  we  have  tested 
this  assumption  in  one  of  the  inversions  below. 

The  inversion  proceeds  using  the  method  described  by  Leeds 
et  al . (1974) . We  calculate  the  phase  travel  time  for  the  n^ 

path  and  the  p^^  period  as 


t V £ . s . 

np  = in  ip 


where  is  the  path  length  of  the  n^  path  in  the  i^  region 

(i=l,...,6)  and  s.  is  the  (theoretical)  phase  slowness  for  the  i^ 
region  at  the  p period.  The  phase  slownesses  s^^  are  functions 
of  the  model  parameters  in  each  region. 


T V (t  - t )• 
n p ""P"  ""P 


where  t is  the  observed  travel-time  for  the  n^^  path  at  the 

npo  ^ ^ 

period.  The  minimization  takes  place  with  respect  to  the  choice 
of  model  parameters. 
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We  obtain  larger  and  larger  variances  in  the  model  para- 
meters as  the  number  of  model  parameters  increases.  We  would 
like  to  be  able  to  solve  for  the  properties  of  the  crust  in  each 
of  the  six  regions.  However,  this  a)  requires  much  precise  data 
at  periods  shorter  than  30  sec.,  b)  increases  the  number  of  model 
parameters  significantly,  and  c)  pushes  our  postulate  of  lateral 
homogeneity  in  each  of  the  regions  to  an  untenable  extreme.  We 
have  therefore  used  model  crusts  for  each  of  the  regions  which  a) 
seem  to  be  plausible  when  compared  with  results  for  similar  parts 
of  the  earth  where  observations  exist  (such  as  locations  typical  of 
regions  1,  2 and  3) , and  b)  agree  with  Soviet  refraction  results 
where  available  (Kosminskaya , et  al,  1969;  Sollogub,  1969).  When 
large  residuals  were  encountered  at  short  periods,  such  as  in  the 
case  of  region  3 (Alpide-Himalayan  belt)  and  region  4 (Tibetan  plateau) , 
we  were  obliged  to  introduce  more  low-velocity  material  into  the  crust. 
This  was  done  by  keeping  crustal  velocities  fixed  and  increasing 
crustal  thickness.  In  these  two  cases,  this  procedure  leads  to 
extraordinarily  thick  crusts.  It  should  be  realized  that  these 
model  crustal  thicknesses  are  consequences  of  the  procedures 
used;  if  we  had  chosen  to  lower  the  crustal  velocities,  the 
thicknesses  would  have  been  less.  We  have  used  a crustal  thick- 
ness of  70  km  in  the  Tibetan  Plateau  (region  4)  a value  not  in- 
consistent with  other  recent  estimates  (Chun  and  Yoshii,  1977). 

The  crustal  models  we  have  used  are  listed  in  Table  5 and  are 


assumed  to  be  fixed  in  the  inversions. 


imdlL  -j 


Model  parameters  for  inversion  of  phase  velocity  data 


CRUST 

LID 

CHANNEL 

SUBCHANNEL 


Thickness 


Depth 
- 0 


B(km/sec)  a(km/sec) 


(different  crustal  models  for  each  region) 
(fixed) 


VAR 


VAR 


94 


94 


63 


138 


139 


420( fixed) 


514 


608 


691 


809 


948 


104 


4.65  8.17 


VAR  7.80 


4.80  8.80 


5.128  9.609 


5.283  9.781 


5.419  9.902 


6.172  10.552 


6.266  11.238 


6.351  11.392 


p(gm/cm^) 

3.45 

3.45 

3.65 

3.806 

3.934 

4.051 

4.417 

4.505 

4.579 


1052 


Crustal  Models  Used  in  Inversion 


crustal  thickness  = 70  crustal  thickness  = crustal  thickness 
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We  have  parameterized  the  mantle  into  a lid,  channel  and 
subchannel;  each  of  these  layers  in  a given  region  is  taken  to  be  homo- 
geneous. The  subchannel  region  terminates  at  a depth  of  420  km. 

Below  this  depth,  we  place  a standard  lower  mantle  platform  under 
all  regions,  indicated  in  Table  5.  The  unknowns  in  the  inversions 
will  be  the  shear  velocities  in  each  of  the  three  upper  mantle 
layers  and  the  depths  of  the  interfaces  at  top  and  bottom  of  the 
second  ("channel")  layer.  The  number  of  unknowns  is  thus  30, 
five  in  each  of  the  six  regions.  This  figure  far  exceeds  the  number 
of  degrees  of  freedom  in  the  data.  Thus  some  reduction  in  the 
number  of  unknowns  has  to  be  made,  by  some  assumptions  which  are 
geophysical  in  character. 

Inversion  1 

In  a first  attempt  at  inversion,  we  have  fixed  the  lid  S-wave 
velocity  at  4.65  km/sec  and  the  subchannel  S-wave  velocity  at  4.8 
km/sec  in  order  to  reduce  the  number  of  degrees  of  freedom.  The 
value  of  4.65  km/sec  for  the  lid  arises  frequently  in  inversions 
for  other  parts  of  the  world.  For  one  case  in  which  a 4.65  km/sec 
lid  was  not  observed,  namely  for  the  Western  United  States  (Biswas 
and  Knopoff,  1974)  in  which  the  subcrustal  material  has  S-wave 
velocity  around  4.3  km/sec,  we  were  able  to  assume  that  a model  with  a 
4.65  km/sec  lid  with  zero  lid  thickness  was  accepted  by  the  inversion 
and  that  the  4.3  km/sec  value  represents  channel  material  rising 
almost  to  the  base  of  the  crust.  Should  the  lid  velocity  in  some 
region  be  less  than  4.65  km/sec  in  reality,  then  crustal  thick- 
nesses can  be  reduced.  A similar  comment  can  be  made  about  the 
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subchannel  velocity:  if  the  subchannel  velocity  should  turn  out  to 
be  less  than  4.8  km/sec  in  some  regions,  channel  thicknesses  will 
be  reduced. 

This  parameterization  thus  includes  only  two  adjustable  model 
parameters  for  each  region.  These  are  the  lid  thickness  and  the 
channel  S-wave  velocity.  A thirteenth  parameter  is  the  sub- 
channel thickness,  which  is  presumed  to  be  uniform  across  the 
entire  continent  and  hence  has  the  same  value  under  each  region. 
Since  the  crustal  thickness  and  the  depth  to  the  420  km  interface 
are  fixed,  the  parameterization  of  lid  and  subchannel  thicknesses 
is  equivalent  to  a parameterization  of  the  depth  below  the  sur- 
face of  the  top  and  bottom  of  the  channel.  This  parameterization 
has  13  degrees  of  freedom. 

After  adjustment  of  the  crust  by  the  procedures  described 
above  (with  interpretation  of  crustal  parameters  according  to 
the  remarks  above) , the  parameterization  and  cross  sections  used 
in  a linearized  inversion  procedure  are  given  in  Table  5.  The 
superficial  sedimentary  layer  that  is  introduced  in  the  crusts 
of  regions  3 and  4 is  designed  to  reduce  the  residuals  at  the 
shortest  periods. 
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The  starting  values  of  the  thirteen  parameters  in  this  linearized 
inverse  are (in  the  usual  units) : 

Table  6 


Region 

^CH 

^LID 

^SUB 

Depth  to 
top  of 
channel 

Depth  to 
bottom  of 
channel 

1 

■ 4.51 

110 

y 

' ; 

155 

f\ 

2 

4.39 

113 

1 

1 

1 

1 

158 

3 

4.30 

54 

j ( 

150 

119 

270 

4 

4.29 

90 

] 

1 

167 

i 

5 

4.08 

92 

j 

137 

6 

4.38 

103 

\ 

/ 

148 

j 

The  thirteen  parameters  in  the  rectangular  box  in  Table  6 are  those  var- 
ied in  the  inversion.  The  standard  errors  of  the  travel-times  in  the 
inversion  were  taken  to  be  the  same  as  those  used  by  Leeds  et  al. 

(1974),  namely  a = Max  (7,  O.IT)  sec.  These  estimates  were  later 
found  to  be  too  large  and  required  modification 

In  the  inversion,  an  iteration  process  has  been  used  in  which 
the  matrix  of  partial  derivatives  G was  recalculated  whenever  we 
moved  into  a new  portion  of  parameter  snace.  This  process  was  con- 
tinued until  we  obtained  convergence  of  the  variable  model  parameters. 

The  thirteen  eigenvalues  of  the  product  matrix  of  partial  derivatives 
T 

G G,  which  were  obtained  in  the  final  stage  of  the  iteration  process, 
are  given  in  Table  7.  Each  eigenvalue  corresponds  to  an  eigenvector 
which, in  every  case  (except  nos.  5 and  6) ,points  in  a direction  close 
to  one  of  the  thirteen  parametric  axes.  Thus,  each  eigenvector  can 
be  said  to  be  a discriminant  for  each  of  the  thirteen  degrees  of 
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TABLE  7 


Eigenvalue 


Model  parameter  most  closely  resolved 
by  eigenvector 


1. 

13.68 

^CH2 

2. 

8.43 

^CH5 

3. 

6.85 

^CH3 

4. 

3.56 

^CHl 

5. 

2.87 

^CH6  ^ ^SUB 

6. 

2.59 

^CH6  ^ ^SUB 

7. 

1.84 

^CH4 

8. 

1.61 

^LID5 

9. 

1.14 

^LID2 

10. 

1.04 

^LID3 

11. 

0.56 

^LID6 

12. 

0.36 

^LIDl 

13. 

0.21 

^LID4 

d 
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> freedom  in  the  model.  The  eigenvectors  corresponding  to  the 

^ largest  eigenvalues  give  reliable  structural  information  con- 

cerning the  model  parameters  closest  to  them,  and  the  eigen- 
vectors corresponding  to  the  smallest  eigenvalues  give  little 
information  about  their  corresponding  parameters. 


i 
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From  Table  7 it  can  be  seen  that,  in  general,  the  data  yield 
the  greatest  information  about  the  channel  velocities  and  the 
least  information  about  the  lid  thicknesses.  In  general,  the 
data  contain  more  information  about  regions  2,  3,  and  5 than 
the  the  other  three  regions.  This  is  clearly  correlated  with 
the  higher  percentages  of  total  path  length  which  sample  regions 
2,  3,  and  5.  (For  complete  details,  see  Table  4 and  Fig.  26). 

The  variance  of  each  model  parameter  decreases  as  successive 
models  approach  the  fine  one,  which  demonstrates  the  convergence 
of  our  inversion  procedure.  Fig.  27  gives  our  results  for  this 
inversion:  the  "best"  model  and  the  corresponding  standard 

deviations  of  each  model  parameter.  We  see  that  the  upper  mantle 
structure  for  regions  i and  2 are  very  similar  although  the  lid 
in  region  1 is  somewhat  thicker.  The  channel  shear-wave  velocity 
for  region  1 is  slightly  less  than  that  in  region  2 , but  the  un- 
certainty of  the  model  parameter  in  region  1 is  rather  large  and 
in  any  case,  the  velocity  contrast  to  the  lid  in  these  regions 
is  rather  smaller  than  one  would  like  to  assert  that  a)  a channel 
is  indeed  present  and  b)  that  a channel  with  velocities  low 
enough  to  require  partial  melting  is  present.  Thus,  the  existence 
of  a low-velocity  channel  in  region  1 is  uncertain.  The  most  strik- 
ing result  of  this  inversion  is  the  very  thin  lid,  and  moderate 
shear  wave  velocity,  in  the  channel  for  regions  3 and  5;  as 
in  the  case  of  regions  1 and  2,  regions  3 and  5 have  strikingly 
similar  upper  mantles.  These  latter  two  regions  are  those  which 
are  tectonically  active  and  have  high  seismicity.  Region  4 
is  that  in  which  the  Eurasian  and  Indian  continents  are 
in  collision.  Although  the  upper  mantle  structure  for  this 
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region  appears  to  be  very  similar  to  those  of  regions  1 and  2, 
the  uncertainties  in  the  model  parameters  are  large.  This  is  due 

1)  to  the  low  percentage  of  total  path  length  in  this  region,  and 

2)  to  the  fact  that  most  of  the  paths  which  sample  this  region 
have  only  short-period  information.  Region  6 has  a thick  lid  and 
a pronounced  low-velocity  channel,  but  the  standard  deviations  are 
rather  large;  additional  assumptions  may  help  to  improve  the 
resolution  of  the  structural  parameters  in  this  region. 

Inversion  2 

In  a second  attempt  at  inversion  we  have  made  the  assumption 
that  the  chemical  composition  of  the  lids  is  everywhere  the 
same  across  Eurasia  (as  before)and  that,  in  this  case,  the  chemical 
composition  of  the  channels  is  similarly  the  same  across  Eurasia.  j 

We  have  allowed  the  lid  velocity  to  be  adjusted  in  this  case;  it  i 

was  fixed  in  the  preceding  case.  In  this  case  the  number  of  para-  j 

meters  in  the  fit  is  nine:  the  depth  to  the  lid-channel  interface  ; 

i 

in  each  of  the  six  regions,  the  depth  to  the  channel  - subchannel 

1 

interface,  and  the  S— wave  velocities  in  the  lid  and  channel;  j 

the  latter  three  parameters  are  are  constants  across  the  j 

entire  region.  The  errors  o are  assumed  to  be  the  same  as  i 

'j 

before.  : 

The  results  of  this  inversion  are  shown  schematically  in  Fig.  28. 

The  result  for  lid  velocity  was  4.56+. 01  km/sec  and  for  channel 
velocity  4.34+. 02  km/sec.  The  results  for  lid  thicknesses  show  striking 
similarity  among  regions  1 , 2,  4 on  the  one  hand  and  3,  5,  and  6 on  the  other. 
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FIG. 


FIG. 


FIG. 


26.  Propagation  path  in  each  subregion  of  the 
Eurasian  continent.  The  percentage  of  total 
path  length  in  each  region  is  shown  in  Table  2. 

27.  Final  result  for  linearized  inversion  of  13 
variable  parameters  in  the  model  and  their 
corresponding  standard  deviations. 

28.  Schematic  upper  mantle  cross-sections  obtained 


in  inversion  2. 


KEV 


REGIONS  Fig.  28 


Lid  S-wave  velocity  = 4,56±0.01  km/sec 
Channel  S-wave  velocity  * 4.34±0.02  km/sec 
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This  observation  lends  support  to  the  interpretation  of  the  first 
inversion,  that  the  sediment-covered  platforms  of  Russia  and  Siberia 
are  largely  sediment-covered  shields.  The  result  for  Tibet  we 
have  obtained  is  that  the  upper  mantle,  to  great  depth,  is 
largely  shield-like  in  character.  We  interpret  this  to  mean, 
if  it  should  be  substantiated,  that  the  Indian  Shield  has  been 
emplaced  beneath  the  Tibetan  crust  during  the  collision  of  the 
plates  and  that  the  mantle  beneath  Tibet  is  presently  at 
temperatures  well  below  the  melting  point  to  great  depth,  i.e. 
the  mantle  under  the  Tibetan  plateau  is  relatively  cool  com- 
pared with  tectonically  active  collision  zones  such  as  region  3. 

On  the  other  hand  regions  3,  5,  and  6 all  exhibit  the  develop- 
ment of  a marked  low-velocity  channel  at  shallow  depths  in  the 
mantle,  which  implies  that  a significant  zone  of  partial  melting, 
perhaps  200  km  thick,  is  present.  All  these  three  regions  are 
thus  presumed  to  be  tectonically  active,  despite  the  low  seismicity 
of  region  6 (the  large  recent  earthquakes  in  China  occurred 
in  the  northern  part  of  region  6) . The  presence  of  a small  channel 
at  great  depth  is  not  considered  to  be  significant:  this  may  be 
an  artifact  of  the  inversion,  due  to  improper  choice  of  parameters 
at  shallower  depths.  These  deep  channels  are  probably  absent 
but  we  cannot  be  absolutely  certain. 

Inversion  2.1 


We  have  made  a variation  on  inversion  2 by  enlarging  the 
extent  of  the  Baltic  shield  in  the  preceding  inversion,  as  shown 
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in  Fig.  29  according  to  a hint  in  the  tectonic  map  of 
Khain  and  Muratov  (1969) . As  might  have  been  expected,  the 
results  of  the  preceding  inversion  are  not  significantly  changed 
in  view  of  the  result  of  inversion  2,  namely  that  regions  1 & 2 
are  virtually  identical  in  cross-section.  Hence  enlarging 
region  1 at  the  expense  of  region  2 cannot  be  expected  to  pro- 
duce a significant  change  in  the  results. 

A tabulation  of  the  results  of  this  inversion  is  as  follows, 
and  is  listed  only  for  purposes  of  comparison  with  the  other 
inversions: 


‘^LID 

= 

4.56+0.01  km/sec 

^CH 

= 

4.34+0.02  km/sec 

^LID 

1' 

= 

206+25  km 

^LID 

2' 

= 

229+29 

^ID 

3 

= 

99+9 

^LID 

4 

= 

246+56 

^LID 

5 

= 

79+7 

^LID 

6 

= 

69+7 

t- 

CH-SUB 

= 

274+8 

Inversion  2.2 

In  each  of  the  two  preceding  inversions  we  have  obtained 
the  result  that  the  upper  mantle  of  both  regions  1 and  2 are 
remarkably  similar.  In  both  cases  we  have  a deep  continental 
root  to  essentially  the  same  depths.  We  have  therefore  made 
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FIG.  29.  Modification  of  regionalization  of  Fig.  25 
by  enlargement  of  region  1 in  Baltic  shield 


area . 
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the  assxunption  that  these  two  regions  are  indistinguishable 
in  a new  inversion  2.2.  We  call  the  coalescence  of  regions 
1 and  2,  a new  region  lA  (Fig.  30) ; there  is  now  no  separate 
region  2. 

As  might  be  expected  the  results  of  the  inversions  are 
not  significantly  changed  from  the  earlier  explorations.  Indeed 
the  results  for  inversion  2.2  are  identical  in  all  respects 
with  that  of  inversion  2.1 


^LID 

= 

4.56+. 01  km/sec 

^CH 

= 

4.34+. 02  km/sec 

^CH-SUB 

= 

274+8  km 

\lD  (lA) 

= 

217+18  km 

\lD{3) 

= 

99+9 

^LID(4) 

= 

246+55 

^LID  (5) 

= 

79+9 

^LID  (6) 

= 

69+7 

all  depths  h are  measured  from  the  surface. 

Inversion  3 

A rather  annoying  aspect  of  the  data  concerns  the  fact  that 
the  phase  travel-time  residuals  from  any  of  the  preceding  models 
are  not  normally  distributed.  There  is  an  unacceptably  large 
number  of  residuals  between  2a  and  4o;  this  result  has  been 
verified  by  a test.  We  have  deleted  several  of  the  data  with 
large  residuals  and  have  proposed  a new  data  set  with  a value  of 
which  places  the  new  data  set  within  acceptable  limits  of  a 
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a normal  distribution  with  a new  value  of  o having  2/3  the 
former  value: 

0=2/3  Max  (7  sec,  O.IT) 

We  have  repeated  inversion  2.2  with  the  revised  data 
set  and  the  new  estimate  of  errors  and  have  obtained  the  follow- 
ing results: 


®LID 

= 

4.57+. 01  km/sec 

^CH 

= 

4.35+. 01  km/sec 

^CH-SUB 

= 

276+5  km 

^LID (lA) 

= 

204+12 

^LID(3) 

= 

99+7 

^LID(4) 

= 

256+41 

^LID (5) 

= 

59+4 

*^LID  (6) 

= 

55+5 

The  shear  velocities  in  the  lid  and  channel  are  unchanged  from  the 
preceding  inversion  2.2.  But  there  have  been  some  significant 
changes  in  the  lid  thickness  of  some  of  the  regions.  Regions  lA 
and  4 continue  to  have  upper  mantles  consistent  with  ancient  (cold) 
shield  models.  However,  regions  5 and  6 now  have  very  thin  lids, 
implying  high  heat  flow  and  the  presence  of  strong  tectonic  pro- 
cesses . 

Inversion  3.1 

To  test  the  validity  of  the  linearized  inversion  procedures 
especially  in  view  of  our  observations  that  the  final  results  are 


strongly  dependent  on  smoothing  procedures  in  the  model  space, 
we  have  applied  the  power  of  our  non-linear  Hedgehog  program 
whose  results  are  independent  of  smoothing  in  the  model  space. 

The  data  set  is  the  reduced  data  set  of  Inversion  3 with  a 
normal  distribution  of  residuals  relative  to  model  2.2.  We 
have  further  fixed  the  channel-subchannel  interface  at  the  value 
given  by  inversion  3.  The  remaining  seven  parameters  were  ex- 
plored in  the  space  given  by  Table  6. 


Parameters  and  Range  of  Search  in  Hedgehog  Inversion 
for  Inversion  3.1 


Starting  Value 

Step  Size 

Lower  Limit 

Upper  Limit 

PI 

4.57 

0.1 

4.47 

4.67 

P2 

4.35 

0.1 

4.25 

4.45 

P3 

204 

30 

144 

264 

P4 

99 

30 

69 

189 

P5 

256 

30 

136 

256 

P6 

59 

30 

59 

179 

P7 

55 

30 

55 

175 

PI: 

Lid  shear  wave  velocity 

in 

Km/Sec 

P2: 

Channel  shear 

wave  velocity  in  Km/Sec 

P3: 

Lid  thickness 

of  region 

lA 

ir 

Km 

P4: 

Lid  thickness 

of  region 

3 

in 

Km 

P5; 

Lid  thickness 

of  region 

4 

in 

Km 

P6: 

Lid  thickness 

of  region 

5 

in 

Km 

P7; 

Lid  thickness 

of  region 

6 

in 

Km 

^ A 
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The  results  of  the  inversion  gave  the  67  acceptable  solutions, 
acceptable  within  1.4o  under  the  postulate  that  there  is  equal 
tradeoff  between  the  effects  of  uncertainties  in  the  model  space 
and  fit  to  the  data  (Table  7) . The  first,  and  most  obvious 
result  is  that  the  lid  velocity  was  accepted  to  be  4,57  km/sec; 
no  other  values  are  acceptable.  The  channel  velocities  accepted 
were  only  4.25  and  4.35  km/sec.  Although  the  lid  thicknesses 
for  the  shields  of  region  lA  could  vary  between  174  and  234  km 
and  for  the  Tibetan  plateau  (region  4)  between  226  and  256  km 
a)  under  some  circumstances  both  lids  could  have  the  same  thick- 
ness (models  14,  30,  31,  etc.)  and  b)  under  no  circumstances 
could  the  Tibetan  plateau  have  a lid  which  was  as  thin  as  (say) 

30  km,  a result  which  would  have  implied  a mantle  appropriate 
to  a tectonically  active  region,  i.e.,  one  with  a well-developed 
high-contrast  low  velocity  channel.  Finally,  although  the 
linear  inverse  of  inversion  3 implied  a difference  in  lid 
thickness  between  region  3 on  the  one  hand  and  regions  6 and  7 
on  the  other,  the  non-linear  model-independent  inverse  gave  some 
solutions  in  which  the  lid  thicknesses  are  such  that  the  lid 
channel  interface  is  roughly  at  a common  depth  below  the 
surfaces  of  all  three  regions  (e.g.  solutions  35,  55,  23,  37, 

46,  57).  We  conclude  from  these  inversion  results  that  cold 
shield  properties  extend  to  great  depth  under  regions  lA  and  4 
while  regions  3,  5,  6 may  have  similar  properties  (at  least 


Table  7 
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Successful  solutions  in  Hedgehog  inversion  3.1 


PI 

P2 

P3 

P4 

P5 

P6 

P7 

1 

4.57 

4.35 

204 

99 

256 

59 

55 

2 

4.57 

4.35 

174 

99 

256 

59 

55 

3 

4.57 

4.35 

234 

99 

256 

59 

55 

k 

4.57 

4.35 

20 /+ 

69 

256 

59 

55 

5 

4.57 

4.35 

204 

1 29 

256 

59 

55 

6 

4.57 

4.35 

204 

99 

226 

59 

55 

7 

4.57 

4.35 

204 

99 

256 

89 

55 

8 

4.57 

4.35 

204 

99 

256 

59 

85 

9 

4.57 

4.35 

174 

69 

256 

59 

55 

10 

4.57 

4.35 

174 

129 

256 

59 

55 

1 1 

4.57 

4.35 

234 

69 

256 

59 

55 

12 

4.57 

4.35 

234 

129 

256 

5^ 

55 

13 

4.57 

4.35 

174 

99 

226 

59 

55 

14 

4.57 

4.35 

234 

99 

226 

59 

55 

15 

4.57 

4.35 

174 

99 

256 

89 

55 

16 

4.57 

4.35 

234 

99 

256 

89 

55 

17 

4.57 

4.35 

174 

99 

256 

59 

35 

18 

4.57 

4. 35 

234 

99 

256 

59 

85 

19 

4.57 

4.35 

204 

69 

226 

59 

55 

20 

4.57 

4.35 

204 

129 

226 

59 

55 

21 

4.57 

4.35 

204 

69 

256 

89 

55 

22 

4.57 

4.35 

204 

129 

256 

89 

55 

23 

4.57 

4.35 

204 

69 

256 

59 

85 

24 

4.57 

4. 35 

204 

129 

256 

59 

85 

25 

4.57 

4.35 

204 

99 

226 

89 

55 

26 

4.57 

4.35 

204 

99 

226 

59 

85 

27 

4.57 

4.35 

204 

99 

256 

89 

85 

28 

4.57 

4.35 

174 

69 

226 

59 

55 

29 

4.57 

4.35 

174 

129 

226 

59 

55 

30 

4.57 

4. 35 

234 

69 

226 

59 

55 

31 

4.57 

4.35 

234 

129 

226 

59 

55 

32 

4.57 

4.35 

174 

69 

256 

89 

55 

33 

4.57 

4.35 

174 

129 

256 

89 

55 

34 

4.57 

4.35 

234 

69 

256 

89 

55 

Lk J 


Inversion  4 


> To  test  the  importance  of  diffraction,  we  have  deleted 

from  the  data  set, phase  travel  times  for  those  paths  which  lie 
, close  and  parallel  to  a contrast  between  two  dissimilar  zones. 

The  value  of  o continues  to  be  taken  as  in  the  result  of 

a, 

Inversion  3.  The  linear  inverse  using  the  same  regionalization 
and  parameterization  as  inversion  3 for  the  results: 


^LID 

= 

4.58+. 01 

km/sec 

^CH 

= 

4.36+. 02 

km/sec 

^LID(IA) 

= 

181+13 

km  (below 

surface) 

^LID (3) 

= 

101+16 

M 

II 

^LID(4) 

= 

204+35 

II 

II 

\lD(5) 

= 

54+7 

II 

II 

^LID  (6) 

= 

54+7 

II 

II 

^CH-SUB 

= 

268+6 

II 

II 

We  detect  no  significant  changes  from  the  results  of  inversion  3. 
Inversion  5 

Finally,  since  the  Tibetan  structure  appears  to  be 
associated  with  the  collision  of  the  Indian-Asian  plates,  we 
have  incorporated  the  Himalayan  part  of  region  3 into  region  4, 
and  re-analyzed  the  linear  inversion  with  the  data  and  errors 
as  in  inversion  3.  A map  of  the  new  regionalization  is  shown 
in  Figure  31. 
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FIG.  31.  Modification  of  regionalization  of  Fig,  30 

by  enlargement  of  region  4 in  Himalayan  region. 


45 


The  results  of  the  inversion  are 


^LID 

= 

4.58+. 01 

^CH 

= 

4.34+. 01 

^LID(IA) 

= 

194+10 

^LID(3) 

=: 

97+8 

\lD{4) 

= 

207+21 

^LID(5) 

= 

67+4 

^LID (6) 

= 

58+5 

CH-SUB 

= 

273+5 

Again  we  observe  no  major  change  in  the  structure. 

We  conclude  from  all  these  tests  that  the  stable  sediment- 
covered  regions  of  the  USSR  are  probably  stable  preCan±>rian 
shields  covered  by  sediments,  that  the  Tibetan  plateau  is  under- 
lain by  relatively  cold  shield  mantle  material  with  no  major 
low-velocity  zone  at  even  moderate  depths  that  might  be  ex- 
pected of  tectonically  active  zones  or  of  young  stable  regions, 
that  the  mountainous  collision  zone  between  the  Asian  and 
Indian  plates  have  a prominent  low  velocity  zone  at  moderate 
depths,  that  the  Sinkiang-Mongolian  active  seismic  zone  has  a 
similar  structure  and  that  South  Eastern  China  is  also  a 
region  of  tectonic  activity  as  indicated  by  the  similarity  of 
its  uppper  mantle  structure  to  the  other  two  seismically 
active  regions. 
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Arctic  Study.  We  have  measured  fundamental  mode  Rayleigh 
waves  over  a number  of  paths  crossing  the  Arctic  Ocean.  We 
have  used  as  sources  four  earthquakes  whose  focal  parameters 
are : 


4. 

Lena  River 

Aug. 

25, 

1964 

13:47:20.6 

78.2°N 

126.6 

°E 

10. 

Eastern  Aleutians(l ) 

Feb. 

6. 

1965 

01  :40:33.2 

53.2°N 

161.9 

°W 

11. 

Eastern  Aleutians(2) 

Feb. 

6, 

1965 

16:50:23.6 

53.3°N 

161  .8 

°W 

13. 

Western  Aleutians 

Feb. 

6, 

1965 

04:02:53 

52.1°N 

175.7 

°E 

For  each  of  these  earthquakes  we  have  obtained  initial  phases 
either  from  the  fault  plane  solution  or  from  the  radiation  pat- 
tern for  Rayleigh  waves  (as  in  the  case  of  the  Lena  River  dis- 
cussed above).  We  have  obtained  phase  velocities  by  the  single- 
station method  for  seven  paths  crossing  the  Arctic  over  the 
period  range  50  to  208  sec.  The  paths  are  shown  in  Figure  32  . 

It  can  be  seen  that  none  of  these  are  purely  oceanic  paths. 
The  shaded  area  outlines  out  estimate  of  the  boundary  between 
the  continental  shelf  and  the  deep  ocean  basins.  The  fraction 
of  the  geometrical  path  that  each  event  has  in  the  oceanic  part 
is  as  follows: 


Event 

Total  path  length 

Oceanic  length 

Fraction  oceanic 

1 . Lena  - ESK 

4816  km 

1904  km 

.40 

2.  Lena  - KT6 

3386 

1823 

.54 

3.  East  Aleut(l)  - KON 

7483 

4603 

.62 

4.  East  Aleut(l)  - KTG 

5933 

916 

.15 

5.  East  Aleut(2)  - KEV 

6347 

2334 

.37 

6.  East  Aleut(2)  - ESK 

7818 

2714 

.35 

7 . West  Aleut.  - KEV 

6263 

844 

.13 

IWoOZl 
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To  reduce  the  available  data  to  information  regarding 
purely  oceanic  paths,  we  have  decided  to  use  the  phase  veloci- 
ty data  for  the  Lena  River  event  recorded  at  KEV  (see  discus- 
sion above  for  Eurasia)  as  a typical  continental  value  and  to 
subtract  these  values  , for  the  appropriate  path  length  con- 
tribution, from  the  phase  delays  observed  for  the  above  7 
path-events.  Unfortunately,  the  two  events  East  A1  eut. (1 ) - KTG 
and  West  Aleut. -KEV  have  such  small  parts  of  their  total  path 
that  are  oceanic  that  we  are  subtracting  two  numbers  of  com- 
parable size  and  the  result  is  quite  unstable.  The  unreliabil- 
ity of  the  oceanic  phase  delay  results  for  these  two  cases  has  obliged 
us  to  exclude  them  from  our  data  set.  Accordingly,  we  have  in- 
vestigated the  inversions  of  the  phase  velocity  results  for 
the  five  remaining  paths.  The  relevant  data  are  given  in  Table  8: 


j 
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TABLE  8 

(Pure)  Oceanic  Phase  Velocities  (km/sec) 


T( sec ) 

Lena-ESK 

Lena-KTG 

E.Al (1 )-K0N 

E.A1(2)-KEV  E 

.A1 (2)-ESK 

208 

4.64 

(4.83) 

4.62 

4.55 

4.52 

192 

4.52 

(4.61) 

4.51 

4.42 

4.42 

167 

4.38 

4.32 

4.38 

4.22 

4.28 

139 

4.21 

4.14 

4.24 

4.09 

4.15 

119 

4.08 

4.06 

4.12 

4.04 

4.11 

100 

4.01 

3.99 

4.08 

3.96 

4.07 

69 

3.95 

3.90 

4.00 

3.91 

4.00 

50 

3.91 

3.82 

3.98 

3.89 

3.92 

With 

so  few 

data,  we  have 

not  been  able 

to  regionalize 

the 
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the  deep  Arctic;  the  number  of  degrees  of  freedom  in  the  data 
is  too  small.  The  best  we  can  do  is  to  consider  the  deep  Arctic 
as  a single  province  and  investigate  the  consequences  of  invert- 
ing an  "average"  phase  velocity  for  the  region.  The  average 
phase  velocity  is  obtained  from  the  above  table  by  weighting  by 
the  oceanic  path  length  in  each  case.  The  result  is  (omitting 
the  quantities  in  parentheses) : 


T (sec) 

c (km/sec) 

208 

4.59 

192 

4.47 

167 

4.32 

139 

4.18 

119 

4.10 

100 

4.03 

69 

3.97 

50 

3.92 

These  results  can  be  compared  with  those  obtained  for  Pacific 
paths  by  Leeds  (1973)  from  inversion  of  trans-Pacif ic  phase 
velocity  data  by  methods  similar  to  those  described  above  for 
trans-Eurasian  paths.  The  pure-age  phase  velocities  for  the 
Pacific  can  be  derived  from  the  cross-sections  resulting  from 
the  inversions;  these  are  shown  in  Fig.  33  for  Pacific  ages 
0-10  my,  20-40  my,  85-110  my.  The  Arctic  data  points  are  shown 
as  circles.  The  Arctic  cross-section  averages  out  to  about  a 30  my 
Pacific  structure.  According  to  <-he  model  of  Parker  and  Oldenburg 
(1973) , the  lid  thickness  as  a function  of  the  age  is 


49 


r 

» 


jj 

z(t)  = 9.4t  km,  with  t the  age  in  my.  Thus,  assuming  that  the 
lid  and  channel  S-wave  velocities  are  of  the  same  order  as  in 
the  Pacific,  the  age  at  which  the  Arctic  began  to  open  is  cal- 
culated to  be  about  70  my  (before  present) , an  unexpectedly  small 
quantity. 

Theoretical  seismograms.  The  main  thrust  of  our  work  with 
time  series  synthesis  has  been  directed  toward  improving  the 
efficiency  of  existing  computational  techniques.  This  improve- 
ment is  required  to  permit  us  to  extend  the  information  contained 
on  the  theoretical  seismograms  down  through  the  period  range 
covered  by  the  long-period  instruments  of  the  WWSSN.  Although  up 
to  the  present  time  we  have  concentrated  on  laterally  homogeneous 
structures,  in  all  other  respects  our  models  of  the  earth  have  been 
highly  realistic;  approximately  200  layers  are  being  used  to  model 
the  radial  heterogeneity  of  the  crust-mantle  system  of  a spherical 
earth,  and  the  intrinsic  attenuation  is  included. 

A summary  of  the  general  methods  we  have  applied  in  our 
computations  is  given  by  Kausel  and  Schwab  (1973) , and  Knopoff 
et  al  (1974) . An  elaboration  of,  and  certain  justifications 
for  these  procedures  have  recently  been  given  by  Schwab  and 
Kausel  (1976) . A recent  contribution  by  Calcagnile  et  al  (1976) , 
also  completed  under  this  contract,  is  also  pertinent  here  when 
only  the  surface-wave  portion  of  the  thoeretical  seismogram  is 
desired. 


The  initial  development  of  the  algorithm  and  certain  pro- 
gramming improvements,  which  were  carried  out  under  the  present 
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' ' FIG.  33.  Results  of  phase  velocity  dispersion  for  Arctic 

region  (circles) . "Standard"  phase  velocity 
curves  for  Pacific  age  strips  are  shown  for 

(- 

comparison  (Leeds,  1973).  Because  the  Pacific 
spreads  at  a different  rate  than  the  Arctic, 
this  gives  a different  age  for  the  Arctic  (see 
text) . 


Fig.  33 
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contract,  are  contained  in  Nakanishi,  Schwa!>  and  Kausel  (1976), 
and  Nakanishi,  Schwab  and  Knopoff  (1976) . These  manuscripts 
formed  the  necessary  bridge  between  first  generation  and  second 
generation  dispersion  programs.  First  generation  techniques 
are  based  on  full,  detailed  specification  — a structure,  a 
specific  mode,  and  a specific  period  — from  which  a single 
phase  velocity  at  a time  is  sought;  second  generation  dispersion 
computations  begin  with  only  the  structure  and  the  mode  specified, 
and  they  then  compute  all  phase  velocities  down  to  whatever 
minimum  period  is  desired.  As  the  second  generation  software  is 
developed,  computations  for  the  group  velocity,  phase  attenuation, 
amplitude  excitation  function,  and  apparent  initial  phase  are  in- 
corporated into  the  procedure  (but  not  yet  into  a single,  automatic 
routine  combined  with  the  phase  velocity  computations) . 

In  a series  of  five  later  papers  under  this  contract  (Kausel, 
Schwab  and  Mantovani,  1977;  Mantovani  et  al  1976,  1977;  Mantovani , 
1977  a,b)  this  second  generation  software  was  fully  developed  and 
applied  to  the  generation  of  multimode  theoretical  siesmograms 
containing  as  many  as  21  modes;  each  of  these  was  represented 
over  the  entire  period  range  down  to  1 second. 

The  final  stage  of  our  work  on  the  generation  of  complete 
theoretical  seismograms  for  torsional  waves,  was  the  development 
of  a third  generation  program.  This  routine  is  fully  automatic, 
and  requires  only  the  structural  specification  as  input.  The  out- 
put, which  is  obtained  in  a single,  relatively  short  computer  run, 
contains  all  of  the  frequency-domain  information  required 
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to  compute  theoretical  seismograms  for  arbitrary  source  specifi- 
cations. As  desired,  all  body-wave  and  surface-wave  arrivals, 

% 

for  periods  greater  than  ten  seconds,  are  obtained  from  the  90-100 
modes  thus  specified.  Results  from  this  work  were  recently  pre- 
pared and  a preprint  is  appended  (Liao,  Schwab  and  Mantovani,  1977). 
VVe  are  now  in  a position  of  being  able  to  compare  directly,  the 
entire  experimental,  torsional-wave  seismograms  from  the  long- 
period  WWSSN  instruments  with  those  computed  from  theoretically 
specified  sources  and  structures. 

Our  work  on  the  algorithm  and  the  programming  on  the  Rayleigh-, 
or  spheroidal-wave  theoretical  seismograms,  began  with  detailed 
analysis  and  improvements  of  the  basic  direct  method  for  handling 
such  calculations  on  a spherical,  gravitating  earth.  In  its 
original  form,  this  method  was  initially  develoepd  in  the  series 
of  papers  by  Hoskins  (1920) , Pekeris  and  Jarosch  (1958) , and 
Alterman,  Jarosch  and  Pekeris  (1959)  ; some  numerical  details  con- 
cerning such  computations  were  given  by  Bolt  and  Dorman  (1961) , 
and  later,  by  Takeuchi  and  Saito  (1972).  In  the  appended  manu- 
script (Schwab  et  al,  1977)  we  describe:  (1)  the  newly  developed 
simplifications  of  the  usual  algorithm,  which  has  made  it  possible  — 
we  believe  for  the  first  time  — to  develop  a direct  algorithm  for 
group  velocity  computation,  that  is  independent  of  the  usual 
appeal  to  variational  techniques;  (2)  the  numerical  problems  that 
are  associated  with  this  type  of  computation,  in  quite  detailed 
form,  in  relation  to  what  has  appeared  previously  in  the  literature; 


(3)  our  optimization  of  the  improved  algorithm,  and  the 
efficiency  of  the  optimization  relative  to  various  similar 
computations  problems. 

Relative  to  the  efficiency  determination,  the  most  per- 
tinent results  affecting  the  work  on  this  contract  are:  (1)  that 
such  computations  — for  Rayleigh-,  or  spheroidal  waves  on  a 
spherical,  gravitating  earth  — are  about  six  times  more  expensive 
than  the  comparable  torsional- wave  computations,  and  (2)  that 
spheroidal  waves  can  be  handled  on  a non-gravitating  earth  for  only 
about  twice  the  expense  of  torsional  waves.  From  these  results  we 
conclude  that  the  present  optimization  is  still  too  expensive 
to  use  for  the  computation  of  theoretical  seismograms  for 
spheroidal  waves  (down  to  periods  of  only  ten  seconds) , but  that 
the  optimization  technique  will  be  satisfactory  for  this  purpose 
if  a means  can  be  devised  to  approximate  the  removal  of  gravity 
from  the  formulation.  Such  a technique  has  already  been  devised 
(Schwab,  1977),  but  the  numerical  tests  have  only  just  begun. 

The  final,  practical  purpose  of  our  work  under  this  contract  — 
application  of  our  results  to  the  discrimination  problem  — will 
involve  comparison  of  complete  theoretical  and  experimental 
seismograms.  It  is  theorefore  important  that  we  have  as  accurate 
a means  as  possible  of  obtaining  the  instrumental  constants  from 
the  impulse  response  of  the  experimental  record.  These  constants 
then  permit  us  to  include,  with  a minimum  of  error,  the  effect 
of  the  instrumental  response  on  the  theoretical  seismogram.  Our 
improved  scheme  for  inversion  of  the  impulse  response  to  obtain 
the  instrument  parameters,  is  described  in  the  appended  preprint 
by  Mitchel  (1977) . 


Bibliography 


53 


Alterman,  Z.,  H.  Jarosch  and  C.  L.  Pekeris  (1959).  Oscillations 
of  the  Earth,  Proc.  Roy.  Soc.,  A 252 , 80-95. 

Alterman,  Z.  and  D.  Loewenthal  (1972).  Computer  generated 

seismogran.3 , in  Methods  in  Computational  Physics,  volume  12 
(edited  by  B.  A.  Bolt),  Academic  Press,  New  York/25-164. 

Arkangel ' sakaya,  V.  M.  (1960).  Dispersion  of  surface  waves  and 

crustal  structure.  Bull.  (Izv.)  Acad.  Sci.  USSR,  Geophys.  Ser. 
No.  9,  904-927. 

Backus,  G.  and  F.  Gilbert  (1968).  The  resolving  power  of  gross 
Earth  data,  Geophys . J.  Roy . astr . Soc.,  169-205. 

Backus,  G.  and  F.  Gilbert  (1970).  Uniqueness  in  the  inver- 
sion of  inaccurate  gross  Earth  data,  Phil . Trans . Roy. 

Soc.  bond. , Ser.  A,  266,  123-192. 

Ben-Menahem,  A.  and  D.  G.  Harkrider  (1964) . Radiation  patterns 
of  seismic  surface  waves  from  buried  dipolar  point  sources 
in  a flat  stratified  earth,  J.  Geophys.  Res.,  69 , 2605-2620. 

Ben-Menaham,  A.  and  M.  N.  Toksoz  (1963).  Source  mechanism  from 
spectra  of  long-period  seismic  surface  waves.  2.  The 
Kamchatka  earthquake  of  November  4,  1952,  J.  Geophys.  Res., 

68,  5207-5222. 

Biswas,  N.  N.  and  L.  Knopoff  (1970).  Exact  earth-flattening 

calculation  for  Love  waves.  Bull.  Sesimol.  Soc.  Amer.,  60, 
1123-1137. 

Biswas,  N.  N.  and  L.  Knopoff  (1974).  Structure  of  the  upper 
mantle  under  the  United  States  from  the  dispersion  of 
Rayleigh  waves,  Geophys.  J.  R.  astr.  Soc.,  36 , 515-539. 


Bolt,  B.  A.  and  J.  Dorman  (1961) . Phase  and  group  velocities 
of  Rayleigh  waves  in  a spherical,  gravitating  Earth, 


54 


J.  Geophys.  Res.  66 , 2965-2981. 

Brune,  J.  and  J.  Dorman  (1964).  Seismic  waves  and  earth 

structure  in  the  Canadian  shield,  Bull.  Seis.  Soc.  Am., 

§2,  167-210. 

Brune,  J. , J.  Nafe  and  J.  Oliver  (1960).  A simplified  method 
for  the  analysis  and  synthesis  of  dispersed  wave  trains. 

J.  Geophys.  Res.,  65 , 287-305. 

Burridge,  R.  and  L.  Knopoff  (1964) . Body  force  equivalents 

for  seismic  dislocations.  Bull.  Seismol.  Soc.  Amer.,  54 , 
1875-1888. 

Calcagnile,  G. , G.  F.  Panza,  F.  Schwab  and  E.  G.  Kausel  (1976). 

On  the  computation  of  theoretical  seismograms  for  multi- 
mode  surface  waves,  Geophys.  J.  R.  astr.  Soc.  47 , 73-81. 

Chun,  Kin-Kip  and  Toshikatsu  Yoshii  (1977) . Crustal  structure 
of  the  Tibetan  plateau:  A surface  wave  study  by  a Moving 
Window.  Bull.  Seis.  Soc.  Amer.,  67 , 735-750. 

Dorman,  J.  and  M.  Ewing  (1962) . Numerical  inversion  of  surface 
wave  dispersion  data  and  the  crust-mantle  structure  in  the 
New  York-Pennsylvania  area,  J.  Geophys.  Res.,  67 , 5227-5241. 

Frez,  J.  and  F.  Schwab  (1976).  Structural  dependence  of  the 
apparent  initial  phase  of  Rayleigh  waves,  Geophys.  J.  R. 
astr.  Soc.,  44 , 311-331. 

Gilbert,  F.  and  A.  M.  Dziewonski  (1975) . An  application  of 


normal  mode  theory  to  the  retrieval  of  structure  parameters 
and  source  mechanisms  from  seismic  spectra,  Phil.  Trans.  Roy. 
Soc.  London,  Ser.  A,  278,  187-269. 


55 

Harkrider,  D.  G.  (1964) . Surface  waves  in  multilayered  elastic 

media.  Part  I.  Rayleigh  and  Love  waves  from  buried  sources 
in  a multilayered  elastic  half-space,  Bull.  Seismol.  Soc.  Amer. 

627-679. 

Harkrider,  D.  G.  (1970) . Surface  waves  in  multilayered  elastic 
media.  Part  II.  Higher  mode  spectra  and  spectral  ratios 
from  point  sources  in  plane  layered  earth  models.  Bull . 

Seismol.  Soc.  Amer.,  60 , 1937-1987. 

Haskell,  N.  A.  (1953) . The  dispersion  of  surface  waves  on 

multilayered  media.  Bull.  Seismol.  Soc.  Amer.,  43,  17-34. 

Hoskins,  L.  M.  (1920) . The  strain  of  gravitating  sphere  of 

variable  density  and  elasticity.  Trans.  Am.  Math.  Soc., 

2^,  1-43. 

Jackson,  D.  D.  (1972) . Interpretation  of  inaccurate,  insuffi- 
cient and  inconsistent  data,  Geophys . J. , 28 , 97-109. 

James,  D.  E.  and  A.  T.  Linde  (1971) . A source  of  major  error 
in  the  digital  analysis  of  World  Wide  Standard  Station 
seismograms.  Bull.  Seismol.  Soc.  Amer. , 61 , 723-728. 

Kausel,  E.  G. , A.  R.  Leeds  and  L.  Knopoff  (1974).  Variation 
of  Rayleigh  wave  phase  velocities  across  the  Pacific 
Science,  186 , 139-141. 

Kausel,  E.  and  F.  Schwab  (1973) . Contributions  to  Love  wave 


transformation  theory:  Earth-flattening  transformation 
for  Love  waves  from  a point  source  in  a sphere.  Bull. 
Seismol.  Soc.  Amer. , 63,  983-993. 


56 


Kausel,  E.  G.,  F.  Schwab  and  E.  Mantovani  (1977).  Oceanic  Sa 
Geophys.  J.  R.  astr.  Soc. , 50 , 407-440. 

Keilis-Borok,  V.  I.  and  T.  B.  Yanovskaya  (1967) . Inverse  seis- 
mic problems  (structural  review) , Geophys.  J. , 223- 

233. 

Khain,  V.  E.  and  M.  V.  Muratov  (1969) . Crustal  movements  and 
tectonic  structure  of  continents,  in  The  Earth's  Crust 
and  Upper  Mantle , (edited  by  P.  J.  Hart) , Geophysical 
Monograph  13,  American  Geophysical  Union,  Washington,  DC. 

Knopoff,  L.  (1961).  Green's  function  for  eigenvalue  problems 

and  the  inversion  of  Love  wave  dispersion  data,  Geophys.  J. , 

4,  161-173. 

Knopoff,  L.  (1962) . Higher  order  Born  approximation  for  the 

inversion  of  Love  wave  dispersion,  Geophys . J.  , 1_,  149-157. 

Knopoff,  L.  (1964).  A matrix  method  for  elastic  wave  problems. 

Bull.  Seismol.  Soc.  Amer.,  54 , 431-438. 

Knopoff,  L.  (1969) . Phase  and  group  slownesses  in  inhomogeneous 
media,  J.  Geophys.  Res.,  74 , 1701. 

Knopoff,  L.  (1972) . Observation  and  inversion  of  surface 
wave  dispersion,  Tectonophysics , 13 , 497-519. 

Knopoff,  L.  and  D.  D.  Jackson  (1972) . The  analysis  of  undetermined 

and  over determined  systems,  in  The  Dynamic  Response  of  Structures 
G.  Herrmann,  ed.,  Pergamon  Press,  289-305. 


57 


Knopoff,  L.,  S.  Mueller,  and  W.  L.  Pliant  (1966).  Structure 
of  the  crust  and  upper  mantle  in  the  Alps  from  the  phase 
velocity  of  Rayleigh  waves.  Bull  Seismol.  Soc.  Amer.  56 , 
1009-1044. 

Knopoff,  L.  and  J.  W.  Schlue  (1972).  Rayleigh  wave  phase 

velocities  for  the  path  Addis  Abba-Nairobi , Tectonophysics , 
157-163. 

Knopoff,  L.  and  F.  Schwab  (1968) . Apparent  initial  phase  of 

a source  of  Rayleigh  waves,  J.  Geophys.  Res.,  73,  755-760. 

Knopoff,  L. , F.  Schwab  and  E.  Kausel  (1973).  Interpretation  of 
Lg,  Geophys.  J.  R.  astr.  Soc.,  33 , 389-404. 

Knopoff,  L. , F.  Schwab,  K.  Nakanishi  and  F.  Chang  (1974) . 

Evaluation  of  Lg  as  a discriminant  among  different  contin- 
ental crustal  structures,  Geophys.  J.  R.  astr.  Soc.,  39 , 

41-70. 

Kosminskaya,  I.  P.,  N.  A.  Belyaevsky  and  I.  S.  Volvovsky  (1969). 
Explosion  seismology  in  the  USSR,  in  The  Earth’s  Crust  and 
Upper  Mantle,  (edited  by  P.  J.  Hart) , Geophysical  Monograph  13 
American  Geophysical  Union,  Washington,  DC. 

Leeds,  A.  R.  (1973).  Rayleigh  wave  dispersion  in  the  Pacific  Basin 
Ph.D.  thesis.  University  of  California,  Los  Angeles. 

Leeds,  A.  R.  (1975) . Lithospheric  thickness  in  the  Western 
Pacific,  Phys.  Earth  Planet.  Int. , 11 , 61-64. 

Leeds,  A.  R. , L.  Knopoff  and  E.  G.  Kausel  (1974).  Variations 

of  upper  mantle  structure  under  the  Pacific  Ocean,  Science , 
186,  141-143. 


Liao,  H.  M. , F.  Schwab  and  E.  Mantovani  (1977).  Computation 
of  complete  theoretical  seismograms  for  torsional  waves 


58 


(appended) . 

Lubimova,  E.  A.  and  B.  G.  Polyak  (1969) . Heat  flow  map  of 
Eurasia,  in  The  Earth's  Crust  and  Upper  Mantle  (edited 
by  P.  J.  Hart) , Geophysical  Monograph  13,  American  Geo- 
physical Union,  Washington,  DC. 

Mantovani,  E. , F.  Schwab,  H.  Liao  and  L.  Knopoff  (1976). 

Generation  of  complete  theoretical  seismograms  for  SH. 

Part  II,  Geophys . J.  R.  astr.  Soc . , 48 , 531-536. 

Mantovani,  E. , F.  Schwab,  L.  Knopoff  and  H.  Liao  (1977).  Tele- 
seismic  Sn:  A guided  wave  in  the  mantle,  Geophys.  J.  R. 
astr.  Soc.  , 709-726  , 

Mantovani,  E.  (1977a).  Generation  of  complete  theoretical 

seismograms  for  SH.  Part  III,  submitted  to  Geophys.  J.  R. 
astr.  Soc.  (appended)^. 

Mantovani,  E.  (1977b).  Generation  of  complete  theoretical  seis- 
mograms for  SH.  Part  IV,  in  preparation. 

Mitchel,  R.  G.  (1977) . Improved  scheme  of  inversion  of  seis- 
mometer impulse  response  to  obtain  electromagnetic 
instruments  constants  (appended). 

Nakanishi,  K. , F.  Schwab  and  E.  G.  Kausel  (1976).  Inter- 
pretation of  Sa  in  a continental  structure,  Geophys.  J.  R. 
astr.  Soc . , 47 , 211-223. 

Nakanishi,  K. , F.  Schwab  and  L.  Knopoff  (1976).  Generation  of 
complete  theoretical  seismograms  for  SH.  Part  I, 

Geophys.  J.  R.  astr.  Soc.,  48 , 525-530. 


59 


Parker,  C.  L.  and  D.  E.  Oldenburg  (1973) . Thermal  model  of 
ocean  ridges.  Nature  Physical  Science,  242 , 137-139. 

Pekeris,  C.  L.  and  H.  Jarosch  (1958) . The  free  oscillations 
of  the  Earth,  in  Contributions  in  Geophysics,  Volume  1, 
Pergamon  Press,  London,  171-192. 

Popov,  I.  I.  (1960) . Dispersion  of  long-period  Love  waves  in 

the  continental  and  oceanic  crust  along  the  path  Indonesia- 
Crimea , Bull.  (Izv.)  Acad.  Sci . USSR,  Geophys.  Ser.  No.  10, 
970-973. 

Press,  F.  (1968).  Earth  models  obtained  by  Monte  Carlo  inversion, 
J.  Geophys.  Res.,  73 , 5223-5233. 

Press,  F.  (1969).  The  sub-oceanic  mantle.  Science , 165 , 174-176. 

Saito,  M.  (1967) . Excitation  of  free  oscillations  and  surface 
waves  by  a point  source  in  a vertically  heterogeneous 
earth,  J.  Geophys.  Res.,  12,  3689-3699. 

Sato,  Y. , T.  Usami  and  M.  Landisman  (1968).  Theoretical  seis- 
mograms of  torsional  disturbances  excited  at  a focus 
within  a heterogeneous  spherical  earth  — Case  of  a 
Gutenberg-Bullen  earth  model.  Bull.  Seismol.  Soc.  Amer., 
133-170. 

Savarensky,  E.  F.,  G.  N.  Bozhko,  T.  I.  Kukhtikova,  A.  B.  Peshkov, 
I.  I.  Popov,  B.  N.  Sheshkov,  O.  I.  Yurkevich  and 
L.  M.  Yudakova  (1969) . On  the  earth  structure  in  some 
regions  of  the  USSR  from  surface  wave  data.  Pure  and 
Applied  Geophysics,  73,  99-119. 


A 


\ 


60 


Savarensky,  E.  F.  and  A.  B.  Peshkov  (1968) . On  the  use  of 
surface  wave  velocities  in  choosing  a model  of  crustal 
structure,  Akad.  Nauk  SSR  Izv. , Fizika  Zemli  No.  10, 

79-87. 

Savarensky,  E.  F.  and  Sh.S.  Ragimov  (1958).  Determining  the 

velocity  of  Rayleigh  waves  and  direction  to  the  epicenter 
by  three  nearby  stations.  Bull.  (Izv.)  Acad.  Sci.  USSR, 
Geophys.  Ser.  No.  12,  866-869. 

Savarensky,  E.  F.  and  Sh.S.  Ragimov  (1959) . On  the  determi- 
nation of  the  average  thickness  of  the  earth's  crust 
from  Rayleigh  wave  group  velocity  measurements.  Bull . 

(Izv.)  Acad.  Sci.  USSR,  Geophys.  Ser.  No.  9,  969-971. 

Savarensky,  E.  F.  and  B.  N.  Schechkov  (1961) . The  structure 
of  the  earth's  crust  in  Siberia  and  in  the  Far  East 
from  Love  and  Rayleigh  wave  dispersion  data.  Bull . (Izv. ) 

Acad  Sci.  USSR,  Geophys.  Ser.  No.  5,  454-456. 

Savarensky,  E.  F.  and  D.  I.  Sikharulidze  (1959).  Crustal 
thickness  determinations  from  measured  Love  wave  dis- 
persion, Bull.  (Izv.)  Acad.  Sci.  USSR,  Geophys.  Ser.  No.  6. 

Savarensky,  E.  F. , 0.  N.  Solov'eva  and  B.  N.  Shechkov  (1959). 

Love  wave  recording  at  seismological  station  Moscow 
and  crustal  structure.  Bull.  (Izv.)  Acad.  Sci.  USSR, 

Geophys.  Ser.  No.  5. 

Schwab,  F.  (1970).  Surface  wave  dispersion  computations: 

Knopoff's  method.  Bull.  Seismol.  Soc.  Amer.,  60 , 

1491-1520. 

Schwab,  F.  (1977).  Surface-wave  dispersion  computations:  Improved 
algorithm  for  Rayleigh  waves  on  a spherical  gravitating  earth » 


in  preparation. 


4] 


61 


Schwab,  F.  and  E.  G.  Kausel  (1976).  Long-period  surface  wave 
seismology;  Love  wave  phase  velocity  and  polar  phase 
shift,  Geophys.  J.  R.  astr.  Soc. , 45 , 407-435. 

Schwab,  F.  and  L.  Knopoff  (1970) . Surface  wave  dispersion 
computations.  Bull.  Seismol.  Soc.  Amer.,  60 , 321-344. 

Schwab,  F.  and  L.  Knopoff  (1971) . Surface  waves  on  multi- 
layered anelastic  media.  Bull.  Seismol.  Soc.  Amer.,  61 , 

893-912. 

Schwab,  F.  and  L.  Knopoff  (1972) . Fast  surface  wave  and 

free  mode  computations.  Chapter  3,  in  Methods  in  Com- 
putational Physics,  Volume  11  (edited  by  B.  A.  Bolt) 

Academic  Press,  New  York,  87-180. 

Schwab,  F.  and  L.  Knopoff  (1973) . Love  waves  and  the  torsional 
free  modes  of  a multilayered  anelastic  sphere.  Bull . 

Seismol.  Soc.  Amer.  63 , 1107-1117. 

Schwab,  F.,  G.  F.  Panza,  H.  M.  Liao,  E.  G.  l^usel~and“j . Frez  (1977). 
Surface-wave  dispersion  computation;  Rayleigh  waves  on 
a spherical  gravitating  earth,  submitted  to  Bull.  Seism. 

Soc.  Am.  (appended) 

Shechkov,  B.  N.  (1961).  Structure  of  the  earth's  crust  in 
Eurasia  from  the  dispersion  of  surface  waves.  Bull. 

(Izv.)  Acad.  Sci.  USSR,  Geophys.  Ser.  No.  5,  450-453. 

Shechkov,  B.  N.  (1964).  Seismic  surface  wave  dispersion  and 

Eurasian  crustal  structure.  Bull.  (Izv.)  Acad.  Sci.  USSR, 
Geophys.  Ser.  No.  3,  183-187. 


I 


62 


Shechkov,  B.  N.  (1970) . Surface  wave  group  velocities  on 

Eurasian  paths,  Akad.  Nauk  SSR  Izv.  Fizika  Zemli  No.  8, 
80-87. 

Shechkov,  B.  N.  and  0.  N.  Solov'eva  (1961).  Group  velocities 
of  Rayleigh  waves  traveling  along  a mixed  continental- 
oceanic  path.  Bull.  (Izv.)  Acad.  Sci . USSR,  Geophys.  Ser. 
No.  8,  768-771 

Sikharulidze , D.  I.  and  R.  K.  Makharadze  (1968).  The  problem 
of  using  surface  waves  in  seismic  exploration,  Akad.  Nauk 
Gruzin.  SSR  Soobshch.  52 , 335-360. 

Sollogub,  V.  B.  (1969) . Seismic  crustal  studies  in  South- 
eastern Europe,  in  The  Earth's  Crust  and  Upper  Mantle, 
(edited  by  P.  J.  Hart) , Geophysical  Monograph  13, 

American  Geophysical  Union,  Washington,  DC. 

Sykes,  L.  R.  (1967) . Mechanism  of  earthquakes  and  nature  of 
faulting  on  the  mid-oceanic  ridges,  J.  Geophys.  Res. 

72,  2131-2153. 

Takeuchi,  H.  and  M.  Saito  (1972).  Siesmic  surface  waves,  in 
Methods  in  Computational  Physics,  Volume  11  (edited  by 
B.  A.  Bolt),  Academic  Press,  New  York,  217-295. 

Thomson,  W.  T.  (1950) . Transmission  of  elastic  waves  through 
a stratified  solid  medium,  J.  Appl.  Phys.,  21,  89-93. 

Wiggins,  R.  A.  (1972).  The  general  linear  inverse  problem; 
Implication  of  surface  waves  and  free  oscillations  for 
earth  structure.  Revs.  Geophys.  and  Space  Phys.,  10 , 
251-285. 


63 


I 

I Zakharova,  A.  I.,  L.  M.  Matasova  and  O.  V.  Soboleva  (1971). 

^ Mechanism  of  the  focus  of  the  main  shock  from  instrumental 

f data,  in  Tashkent  Earthquake  of  26  April  1966  (ed.  G.  A. 

I Mavlyanov)  F.A.N.  Tashkent  Publ.  House,  Tashkent,  USSR, 

53-58. 


SURFACE-WAVE  DISPERSION  COMPUTATIONS: 


RAYLEIGH  WAVES  ON  A SPHERICAL  GRAVITATING  EARTH 

1 2 

F.  Schwab,  G.  F.  Panza  , H.  M.  Liao,  E.  G.  Kausel  and  J.  Frez 


Publication  No.  1700 

Institute  of  Geophysics  and  Planetary  Physics 
University  of  California,  Los  Angeles 

Publication  No.  121 
Department  of  Geophysics 

University  of  Chile,  Santiago 


Sponsored  by 

Advanced  Research  Projects  Agency  (DOD) 

ARPA  Order  No.  3291 

Monitored  by  AFOSR  Under  Contract  //F49620-76-C-0038 

^Istitute  di  Geodesia  e Geofisica,  Universita  di  Bari,  Italia; 
Dipartimento  di  Scienze  della  Terra,  Universita  degli  Studi 
della  Calabria,  87030  Castigloni  Scalo--Cas.  Post.  14,  Cosenza, 
Italia . 

2 

Departamento  de  Geofisica,  Universidad  de  Cliile,  Casilla  2777  , 


and 


Santiago,  Chile. 


ABSTRACT 


Algorithmic  and  numerical  analyses  are  carried  out 
for  Ray lelgh-wave  dispersion  computations  on  a spherical, 
gravitating  earth.  Our  work  Is  based  on  the  direct,  Alterman- 
Jar osch-Pekerls  formulation.  For  practical  purposes,  we  fix 
period  and  determine  the  associated  phase  velocity  (or  polar 
order  number).  Neither  this,  nor  integration  downward  from  the 
free  surface  — both  "non-standard"  procedures  — results  in 
unexpected  difficulties.  The  latter  procedure  yields  a simpli- 
fication of  the  computational  algorithm,  the  clarity  of  which 
allows  it  to  be  extended  to  group-velocity  evaulations.  The  AJP, 
direct-integration  formulation  is  optimized  and  compared  with  the 
fastest  --  Knopoff's  method  --  of  the  techniques  based  on  the  flat, 
homogeneous-layer  approximation.  The  optimized  form  of  the  AJP 
method  (spherical)  is  three  times  slower  than  Knopoff's  (flat, 
non-gravitating)  method  when  gravity  is  included  in  the  AJP 
formulation;  and  is  1.36  times  slower  when  gravity  is  not  included. 
Additional  programming  would  reduce  the  former  estimate  to  a lower 
bound  of  2.42  times  slower,  and  the  latter,  to  a lower  bound  of  1.30. 
In  size  and  number,  the  treatment  of  integration  "steps"  in  the 
direct-integration  procedure,  is  equivalent  to  the  treatment  of 
"layers"  in  the  homogeneous-layer  approximation;  thus  the  usual 
assumption  that  the  former  method  does  a better  job  of  treating 
continuous  parameter-depth  distributions,  appears  to  be  invalid. 
Overflow  problems  in  the  AJP  formulation  can  be  controlled  by 


simple  normalization.  Loss-of-precision  problems  appear  to  be 
intrinsic  to  the  AJP  formulation.  At  a fixed  period,  this  results 
in  the  attainable  accuracy  of  the  phase  velocity  decreasing  as 
mode  number  Increases;  and,  for  fixed  accuracy  in  the  phase 
velocity,  as  period  decreases  the  maximum  mode  number  that 
can  be  treated  successfully  decreases. 


1. 


INTRODUCTION 


Dorman,  Ewing  and  Oliver  (1960)  described  the  use  of  an  elec- 
tronic computer  to  calculate  surface-wave  dispersion  for  multilayered 
pe r f ec t ly-elas t ic  half-spaces.  Their  computations  were  based  on 
the  technique  devised  by  Thomson  (1950)  and  Haskell  (1953).  Press, 
Harkrider  and  Seafeldt  (1961)  also  used  the  Thomson-Haskell  tech- 
nique, and  with  a more  advanced  computer,  greatly  improved  the 
speed  of  computation.  Randall  (1967)  later  applied  Knopoff's 
(1964)  method  to  this  problem  and  reported  a further  improvement 
in  speed  for  the  Rayleigh-wave  case. 

In  a later  series  of  papers,  Schwab  (1970)  and  Schwab  and 
Knopoff  (1970,  1971,  1972,  1973)  improved  the  optimization,  for 
computer  application,  of  both  the  Thomson-Haskell  and  Knopoff's 
methods  for  flat,  multilayered  media.  These  papers  also  provide 
complete  details  for  obtaining  full  control  over  the  accuracy  of 
the  computations,  and  for  generalizing  the  algorithms  to  include 
computation  of  attenuation  due  to  the  intrinsic  anelasticity  of 
the  earth . 

For  Love  waves,  the  use  of  spherical-to-f lat  structure  trans- 
formations (Biswas  and  Knopoff,  1970;  Schwab  and  Knopoff,  1971; 

1972;  1973;  Kausel  and  Schwab,  1973)  makes  it  possible  to  carry 
out  all  spherical  dispersion,  attenuation,  and  focal-response 
problems  using  the  optimized  algorithms  for  flat  structures.  Sev- 
eral attempts  have  been  made  to  develope  the  same  type  of  transfor- 
mation for  Rayleigh-wave  computations  (Alterman,  Jarosch  and 
Pekeris,  1961;  Bolt  and  Dorman,  1961;  Biswas,  1972;  Schwab  and 
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and  Knopoff,  1972),  but  these  have  all  yielded  only  empirical  re- 
sults which  lack  general  applicability.  Thus,  at  the  present  time 
at  least,  it  appears  that  one  cannot  apply  transformation  theory  to 
Rayleigh-wave  dispersion  co..iputations  on  any  arbitrary,  spherical, 
gravitating  earth.  Bhat tacharya ' s ( 19 76)  recent  results  --  although 
we  will  not  pursue  this  approach  in  the  present  paper  --  suggest  the 
feasibility  of  an  interesting  new  procedure  for  treating  spherical, 
gravitating  structures:  Gravitation  alone  might  be  handled  by  trans 
formation  techniques,  while  Bhat tacharya ' s approach  could  be  used  to 
optimize  the  treatment  of  sphericity. 

The  primary  purposes  of  the  present  paper  are:  (1)  to  report 

on  our  study  of  the  optimization  of  the  direct  computations  (see 
Wiggins  (1976)  for  a discussion  of  computations  based  on  the 
variational  technique),  (2)  to  report  the  results  of  our  study  con- 
cerning accuracy  considerations,  and  (3)  to  determine  the  efficiency 
of  these  direct  computations  relative  to  the  analogous  computations 
for  non-gravitating  structures.  Also,  a new  computational  technique 
is  developed,  for  the  calculation  of  group  velocities,  which  does 
not  depend  on  the  numerical  evaluation  of  "energy  integrals."  Our 
second  purpose  is  to  present  — we  believe  for  the  first  time  — an 
explicit,  quantitative  comparison  of  the  relative  efficiencies  of 
the  two  basic  techniques  for  performing  surface-wave  dispersion 
computations:  that  in  which  an  exact  structural  specification  is 

employed  with  approximate  mathematical  methods,  and  that  in  which 
exact  analytical  techniques  are  applied  to  an  approximate  model  of 
the  structure,  l.e.,  where  the  structure  is  replaced  by  a sequence 
of  homogeneous  layers. 


3. 


2.  ALTERMAN-JAROSCH-PEKERIS  (AJP)  FORMULATION 


The  basic  formulation  for  our  problem  (Pekeris 
and  Jarosch,  1958)  involves  the  solution  of  three  second-order, 
ordinary  differential  equations  constrained  by  a set  of  boundary 
conditions.  For  purposes  of  numerical  solution  it  is  advisable 
to  reduce  this  system  to  six,  linear,  first-order  differential 
equations,  as  was  done  by  Alterman,  Jarosch  and  Pekeris  (1959). 

Bolt  and  Dorman  (1961)  applied  this  formulation,  to  the  evaluation 
of  Rayleigh-wave  dispersion,  and  reported  on  those  numerical 
details  which  it  was  economically  feasible  to  investigate  with 
second-generation  computing  equipment.  Detailed  algorithmic  testing 
of  accuracy,  precision,  and  efficiency  characteristics  really 
requires  the  present,  third-generation  machinery,  which  we  have 
employed  in  the  current  study;  the  work  we  report  here  can  be 
considered  as  the  logical  extension  of  the  above  series  of  papers. 


To  sketch  the  AJP  formulation,  if  we  let 
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with  Yi  and  related  to  the  components  of  displace- 
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For  a treatment  of  the  situations  which  require 
the  use  of  (2.03),  (2.04),  or  their  sum,  see  Schwab  and  Kausel 

(1976).  In  this  same  reference,  the  justification  is  given  for  our 
major  departure  from  previously  reported  computations  of  Rayleigh 
wave  dispersion  on  a sphere:  Strictly  speaking,  Rayleigh  waves 
only  exist  on  a sphere  at  the  discrete  set  of  frequencies  corres- 
ponding to  integral  values  of  the  polar  order  number  ^ . However, 
fixing  and  evaluating  the  corresponding  angular  frequency  m does 
not  yield  the  dispersion  data  at  equal  frequency  intervals,  which 
we  desire  to  use  in  the  usual  numerical  technique  for  obtaining 
time  series  by  inverse  Fourier  transformation.  Schwab  and  Kausel 
(1976)  have  shown  that,  for  most  practical  applications  of  propaga- 
ting surface  waves,  non-integral  ^ at  equally-spaced  frequencies 
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5. 


can  be  used  without  introducing  significant  errors;  therefore, 
we  adopt  the  procedure  of  fixing  ^ and  computing  or 

c = au)/(£+l/2)  , (2.05) 

where  a is  the  radius  of  the  earth.  The  relation  between  £ 
and  the  true  spherical  phase  velocity  is  also  treated  by  Schwab 
and  Kausel  (1976).  In  equation  (2.01),y2  and  yt^  are,  respectively, 

the  radial  dependences  of  the  £r,  and  the  r^  and  r£  components  of 
stress;  y^  and  yg  arise  from  the  presence  of  the  gravitational 

field  . 


i 

[ 

r 
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Since,  in  any  numerical  integration  procedure, 
it  is  important  to  initiate  the  integration  with  accurate  values, 
we  have  chosen  to  proceed  from  the  free  surface  downward.  This 
allows  us  to  specify  the  initial  vector  exactly.  The  integration 
is  then  carried  down  to  a depth  sufficient  to  make  it  immaterial-- 

to  n significant  figures  in  ^ or  £ just  how  we  terminate  the 

integration:  with  an  approximation  of  a free  surface  or  rigid 

surface,  for  example.  The  fact  that  such  a termination  process  is 
valid  has  been  checked  by  extensive  numerical  tests  in  the  course 
of  this  work.  These  tests  follow  the  lines  of  the  layer-reduction 
experiments  described  by  Schwab  and  Knopoff  (1970;  1972),  and  will  be 
described  in  some  detail  below.  In  Section  7 we  discuss  the  termina- 
tion of  the  structure  at  depth  by  either  a solid  or  liquid,  homogen- 
eous, gravitating,  sphere. 


J 


6. 


Here,  we  should  point  out  that  the  warning 
given  by  Takeuchi  and  Saito  (page  241,  1972)  against  proceeding 
downward  from  the  free  surface  when  integrating  the  system  of 
differential  equations,  or  when  forming  the  layer-matrix  product 
if  applying  the  Thomson-Haskell  technique  or  Knopoff's  method, 
does  not  appear  to  be  justified  by  our  experience.  In  the  work 
upon  which  we  report  herein,  downward  integration  did  not  give 
rise  to  any  unexpected  difficulties;  in  previous,  extensive  work 
with  matrix  methods  applied  to  Rayleigh-wave  dispersion  computations 
(Schwab,  1970;  Schwab  and  Knopoff,  1970;  1972),  the  formation  of 
matrix  products  upward  toward  the  free  surface  (Thomson-Haskell 
formulation)  was  not  found  to  have  any  advantage  over  formation  of 
the  product  in  the  downward  direction  (Knopoff's  formulation). 


Continental  structure.  In  this  case,y,  and 
2 

y vanish,  and  y -- y ( JM-1 ) / a at  r •»  a.  Thus  we  can  write  the  starting 

— 


vector  as 


or 


Yg(a)=yi(a)Xi(a)+y3(a)X2(a)+y5(a)X3(a),  (2.07) 

and  for 

Ys  (r)=yi (a)Xi (r)+y3 (a)X2 (r)+y5 (a)X3 (r)  . (2.08) 

The  three  quantities  which  are  unknown  --  yi(a),  y3(a), 
and  y5(a)  --  can  be  carried  implicitly  in  the  computations, 

while  we  integrate  the  vectors  whose  starting  values  are  known 
exactly:  Xj , X2 > and  X3.  That  is,  we 

integrate  to  obtain  X^  at  depth;  this  is  repeated,  in  turn, with 
X2  and  X3.  Thus  we  actually  use  e qua t ion ( 2 . 01 ) in  the  form 


X =AX  to  integrate  from  the  surface  r = a to  the  depth  at 
i i ~ 

which  the  boundary  conditions  are  to  be  applied:  r = r^,  where 
we  can  again  express  Yg  in  terms  of  the  undetermined  coefficients 

by  using  equation  (2.08). 


If  we  define  a rigid  boundary  at  depth  by 


yi(to)=y 3(^0) =y 5(^0) =0. 


(2.09) 


we  then  obtain  three  linear,  homogeneous  equations  in  three 

unknowns the  undetermined  coefficients and  the  determinant  of 

the  coefficient  matrix  must  vanish  if  we  are  to  have  a non-trivial 
solution.  Thus  the  dispersion  function,  P^,  takes  the  form 


F^(c,a))  = 


[Xi(ro)] 
[Xi(ro)] 
[Xi (ro)] 


[X2(ro)]  1 
[X2(ro)l3 
[X2 (r  0 )] 5 


[X3 (ro)]  1 
[X3 (ro)] 3 
[x  3 ( 0 ) ] 5 


, (2.10) 


zeros  of  which  define  valid  (c,a))  dispersion  pairs.  For  the  two 
approximations  to  free  boundaries  at  depth,  we  have  used  the 
def initions 


y2  (*^o)  = y4  (^o)  “ 0 

and  either 


(2.11) 


ye^’^o)  — y5(»^o) 


(2.12) 
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or 


yeCr# )=-y5(ro) (^+1) /a  . 

which  yield,  respectively,  dispersion  functions 


(2.13) 
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(2.15) 


Oceanic  Structure.  In  this  case,  the  analog  of  equation 
(2.01)  is,  for  the  homogeneous  oceanic  (liquid)  layer, 
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At  r^  • a,  y2  vanishes  and  ye  =“y  5 (^  + 1 ) /a , and  we  can  write  the 
starting  vector  as 


yi 

y2 

ys 

ye 

(2.16) 
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ys  (a) 


ys  (a) 


yg  (a)  -ys  (a)  (^-  + 1)  /a 


yi  (a) 


-(£+1) /a 


(2.17) 


Y,  (a)  = yi(a)  Zi(a)  + ysCa)  Z2(a), 


and  for  r < a 


Yj^(r)  = yi(a)  ZjCr)  + ysCa)  Z2(r). 


(2.18) 


(2.19) 


Again,  we  carry  the  unknown  quantities  --  yi  (a)  and  y5(a)-- 
implicity,  and  integrate  the  vectors  whose  starting  values  we 
know  exactly:  Zj  and  Z2,  using  equation  (2.16)  in  the  form  Z 

On  the  oceanic  side  of  the  liquid-solid  boundary  at  the  bottom 


of  the  ocean,  ^ = ri,we  then  have 
[Z  1 ( r 1 ) J 1 

Y,  (ri)=yi(a)  , , + y5(a) 

^ [Zi(ri)j2 

[Zi (r  1 ) J 5 
[Z  1 (r  1 )]  g 


Cz  2 ( >^  1 ) I 1 

[ z 2 ( >^  1 ) ] 2 
[Z  2 ( r 1 ) ] 5 
[Z  2 ( r 1 ) ] g 


(2.20) 


At  this  boundary,  yj,  y2,  ys.  and  yg  are  continuous,  must 

vanish,  and  y3  is  undetermined.  Thus  on  the  solid  side  of  this 
interface,  we  have 
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Copy 
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3.  ALGORITHM  FOR  GROUP  VELOCITY  DETERMINATION 

When  treating  Rayleigh  waves  on  a spherical,  gravitating 
earth,  the  variational  technique  is  usually  employed  to  compute 
group  velocities  (Takeuchi  and  Saito,  1972,  Section  III).  This 
involves  the  evaluation,  by  numerical  approximation,  of  integrals 
having  the  form 

•a 

f(r)  y.(r)y^(r)  dr  . (3.01) 

•*0 

Since  the  functions  yj^(r)  become  highly  oscillatory  for  large 
(radial)  mode  numbers,  this  numerical  evaluation  can  become 
inaccurate  (Knopoff  e_t  a^.  , 1974,  Appendix).  Further,! 

become  spuriously  large  at  depths  much  below  those  at  which  there 
is  significant  energy  in  the  mode,  at  the  period  being  treated. 
Although,  somewhat  surprisingly,  these  spurious  magnitudes  do  not 
seem  to  affect  the  location  of  a root  of  the  dispersion  function, 
they  can  cause  large  errors  in  the  evaluation  of  integrals  such  as 
(3.01).  Thus  one  must  specify  r^  , the  value  of  £ below  which 
the  spurious  magnitudes  occur,  ^Fior  to  seeking  the  group  velocity, 

and  then  evaluate 

•a 

f (r) y^ (r)  y (r)  dr  (3.02) 

Jr  2 ^ 

in  place  of  (3.01).  However  without  prior  knowledge  of  the  group 
velocity,  tj  can  be  quite  d ..cult  to  determine  in  period  ranges 
such  as  those  in  which  the  e rgy  shifts  back  and  forth  repeatedly 
between  the  crustal  wave  guide  and  the  low-velocity  channel  in  the 
upper  mantle  (Panza,  Schwab  and  Knopoff,  1972;  Frantsuzova,  Levshin 
and  Shkad inskaya , 1972;  Schwab  and  Knopoff,  1971;  1972).  Since 

the  integrals  must  eventually  be  evaluated  to  obtain  excitation 
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functions  for  earthquake  sources,  and  since  group  velocity  can  be 
used  to  specify  whether  the  energy  is  in  the  crust  or  low-velocity 
channel,  it  is  desirable  to  develope  an  algorithm  for  obtaining  these 
velocities,  which  is  not  dependent  upon  prior  knowledge  of 


Since  we  have  a dispersion  function,  F(c,a)),  which  vanishes 
when  the  point  ( c , la ) falls  on  a dispersion  curve,  we  can  use  implicit- 
function  theory  to  define  the  group  velocity  u ; 


u 


c/(l- 


0)  dc^ 
c (i) 


with 


(3.03) 


(3.04) 


where  the  partial  derivatives  are  evaluated  at  a point  on  a 
dispersion  curve.  If  we  agree  to  use  the  rigid  boundary  at  depth, 
then 


[xi(ro)]i[xi(ro)]i  [x^(ro)]i 

[Xi(ro)]i  [X2(ro)]i  [X3(ro)]i 

’"-(iDc- 

[X](ro)j3[X2(ro)J3  [X3(ro)j3 

[Xi  (ro)]  5 [X2(ro)]  5 [X3(ro)]5 

+ 

[xi(ro)]3  [xHro)]3  [xUro)]3 

[Xi(ro)]5  [X2(ro)]5  [X3(ro)]5 

+ 


[Xi(ro)]i[X2(ro)]i  [X3(ro)]i 
[Xi(ro)]3[X2(ro)]3  [X3(ro)]3 
[x|(ro)]5[X2(ro)]5  [X3(ro)]5 


(3.05) 


4 


* , 3 F 

F - (^) 


(3.06) 


is  given  by  the  same  type  of  expression  as  (3.05),  with  dots 
replacing  the  primes.  The  elements  lx^(ro)lj  in  (3.05)  and 
(3.06)  are  obtained  exactly  as  described  in  Section  2.  The 
evaluation  of  |xl(rQ)|j  and  |x^(ro)|j  requires  a simple  extension  of 

the  algorithm. 


Continental  Structure.  In  this  case  we  start  with  the 


sixth-order  system 


and  form 


Y = AY  , 


Y '=A' Y+AY ' 


(3.07) 


(3.08) 


Y=  AY+AY.  (3.09) 

Here  again,  we  use  these  equations  of  motion  in  terms  of  the  vectors. 


w 

X'.,  X^,  that  we  know  exactly  at  the  surface: 


jr  . = A'  X.  + AX' . 
1 11 


X.  = A X.  + AX. 
1 11 


(3.10) 


(3.11) 


Since  X can  be  determined  independently,  we  can  treat  A'X.  and 
i ±. 

A X. 

1 


as  known  vectors  at  each  depth,  and  we  have 


X'  . = AX'  . + 
1.  1 


X.  = AX.  + D.  , 
1 11 


where 


C^(r)=A'  (r)X^(r)=-2p(r)6J 


[Xi(r)]i 


[X.(r)]3 


(3.13) 


(3.14) 


D.  (r)=A(r)X.  (r)=-  §^(2S,  + 1) 
1 1 c ^ r 


[x/(X  + 2m)]  [X.(r)]3 

[pg-2ii(3X+2p)/  (X  + 2y)r]  [x.  (r)]  3+[x^(r)] 

0 

[4ii  (X+ii)  / (X+2u)r]  [X^(r)]3 

0 

[-4TrGp]  [X^(r)]  3+[l/r]  [x^(r)]5 


(3.15) 
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Oceanic  structure.  Here  we  begin  with  the  fourth-order 
system  (2.16),  and  form 


— ~ 

“ — 

— — 

7i 

yi 

y'l 

y*2 

= B' 

yi 

+B 

y'2 

^5 

ys 

y's 

^6 

ye 

j^e 

— — 

— — 

4 

^1 

^2 

% 

^2 

^2 

= B 

+B 

• 

^5 

« 

^6 

^6 

• 

I. 

— « 

(3.16) 


(3.17) 


These  equations  are  then  used  with  the  vectors  for  which  we  have 
solutions  at  the  surface,  Z^: 

Z'.  = B'  Z . + BZ'.  (3.18) 

1 11 

Z . = BZ . + BZ  . , (3.19) 

111 
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which  can  be  written  in  terms  of  known  vectors,  B'Z.  and  BZ . , 

3 X 


at  each  depth: 


Z'.  = BZ'.  + E . 
1 1 1. 


Z . = BZ.  + G.  , 
1 11 


where 


-H.  (r) 
1 


I^(r) 


E^(r)=B'  (r)Z^(r)  = [2)l(J!,+  l)  / r2a)3] 


[4TTGp(r)]  H^(r) 


H^(r) 


G^(r)=B(r)Z^(r)  = [-a(2il+l)/c^r^a)3 


[g(r)p (r^H^(r) 


J.  (r) 


H.(r)-g(r)  [Z.(r)]i-[l/p(r)]  [Z.(r)]2  - [z.(r)]5 


(3.20) 


(3.21) 


(3.22) 


(3.23) 


(3.24) 
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[!  * 
li 

I? 


kk. 


I^(r)=-p(r)  ^‘‘r^ /«,(£+!)  +g^(r)]  [z^(r)]  i + g(r)  { [z^(r)]  2 -H) (r)  [z^(r)]  5}  (3.25) 

J^(r)=4TiG{-p(r)g(r)  [z^(r)]i  +[z^(r)1  2+[p (r)-<<i)^/4TrG]  [z^(r)]5}  (3.26) 

Application  of  (3.20)  and  (3.21)  will  allow  us  to  carry  the  integration  to  the 


bottom  of  the  oceanic  layer  at  where  we  can  apply  the  boundary  conditions 


of  continuity  of  yi,  y'  i.yj,  ya.y' 2.y2.y5.y' 5.y5»y6»y' 6 »y6 » ^nd  the  vanishing  of 


y4»y' 4»y4»  to  obtain  the  necessary  starting  values  to  apply  in  (3.12)  and  (3.13): 


X'l  (ri) 


[z'  1 (ri)]  1 

[Z'2(r)]i 

[2'  1(^1)]  2 

[Z'  2(r)]2 

0 

,X'2(ri)=0,X'3  (ri)- 

0 

0 

0 

[Z*  1 (ri)]5 

[Z'2  (r)]5 

[Z’l  (ri)]6 

[Z’2(r)]6 

(3.27) 


and  a like  set  of  starting  values  with  primes  replaced  by  dots. 

Although  precise  computat ion-t ime  estimates  are  given  in 
Section  5,  a few  general  observations  will  be  appropriate  here.  At 
any  given  frequency,  the  phase  velocity  is  evaluated  by  repeated 
determination  of  the  dispersion  function  until  a root  is  bracketed; 
this  root  Is  then  refined  until  the  desired  accuracy  in  £ is  obtained. 
Each  determination  of  involves  the  integration  of  the  three 


vectors  X..  If  we  assume  that  X.(r)  have  been  saved  from  the  last 
1 1 


set  of  integrations  to  obtain  a value  of  c^,  then  to  obtain  the  group 
velocity  at  this  same  frequency  requires  the  integration  of  six 
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vectors,  i.e.  approximately  twice  the  amount  of  time  it  takes  to 

compute  a single  value  of  the  dispersion  function  while  iterating 

for  c_.  The  number  of  iterations  required  in  the  computation  of  c^ 
is  variable,  but  for  the  present  purpose,  ten  iterations  may  be 

taken  as  a representative  value.  Thus  30  integrations  would  be 

required  to  obtain  c_,  and  only  6 to  obtain  To  the  accuracy  of 

this  estimate  then,  ^ can  be  obtained  in  only  20  percent  of  the  time 

the  evaluation  of  c^  requires. 
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4 . NUMERICAL  TECHNIQUE  FOR  INTEGRATING  THE  SYSTEM  OF  DIFFERENTIAL 
EQUATIONS 

To  optimize  the  computations,  it  is  useful  to  constrain 
the  structural  specification  somewhat.  For  this  purpose: 

1.  The  liquid,  oceanic  layer  is  limited  to  a single, 
homogeneous  layer.  A special  fourth-order  Runge-Kutta  technique 
(see  below)  is  used  for  the  first  three  steps  of  the  integration-- 
step  size  of  about  1 km — and  a fourth-order  predicator-corrector  method  (see 
below)  is  employed,  if  necessary,  to  continue  the  integration  to 
the  bottom  of  the  oceanic  layer  with  the  same  step  size. 

2.  The  sedimentary  layers  are  limited  to  a sequence  of 
homogeneous  layers,  each  of  which  does  not  exceed  1 km  in  thickness. 
These  layers  are  treated  with  one  fourth-order  Runge-Kutta  step. 

3-  The  sub s ed imen t ar y crustal  layers  must  also 
be  homogeneous,  and  each  is  treated  as  the  oceanic  layer;  the  step 
size  fixed  at  about  1 km. 

4.  The  sub-Moho  mantle  requires  continuous  velocity- 
depth  and  density-depth  distributions  , although  discontinuities 


± 


— ■ 


r 


can  be  approximated  as  closely  as  desired  by  specifying  large 
gradients.  As  with  the  oceanic  layer, three  Runge-Kutta  steps 
are  followed  by  a fourth-order  predictor-corrector  method.  The 
initial  step  size  is  1.5625  km,  with  which  we  execute  three  Runge- 
Kutta,  and  seven  predictor-corrector  steps.  The  step  size  is  then 
doubled  and  five  predictor-corrector  steps  are  performed;  this 
procedure  is  repeated  until  the  step  size  reaches  12.5  km,  and  the 
predictor-corrector  method  is  then  applied  with  this  fixed  step 
size.  The  results  of  our  numerical  testing  have  shown  that  this 
dependence  of  step  size  on  depth  is  sufficient  to  yield  4-signif^ 
cant-figure  accuracy  in  the  computed  values  of  c^.  Concerning  this 
point,  one  should  review  the  treatment  given  by  Schwab  and  Knopoff 
(1972),  in  which  piecewise-cont inuous  velocity-  and  density-depth 
distributions  are  treated  with  the  homogeneous-layer  approximation. 
Comparison  will  show  that  the  thickness  of  the  layers  as  a function 
of  depth,  in  that  approximation,  is  roughly  the  same  as  the  integration 
step  size  in  the  present  analysis.  That  is,  to  the  degree  of  accuracy 
which  two  such  dissimilar  methods  can  be  compared,  if  the  step  sizes  above 
are  used  as  layer  thicknesses  in  the  homogeneous  layer  approximation,  that 

technique  will  yield  4 significant  figures  in  the  computed  values  of  c_  . 

f- 

> 

t 
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Runge-Kutta  technique.  To  start  the  predictor-corrector 
method  we  have  used  a Runge-Kutta  technique  designed  specifically 
for  this  purpose.  Here  one  is  only  interested  in  being  able  to 
minimize  the  bounds  on  the  truncation  error.  Ralston  (1962)  has 
treated  this  problem,  and  gives  the  following  algorithm  for 
obtaining  the  first  four  points  for  starting  a predictor- 
corrector  method. 

In  terms  of  a single,  first-order  differential  equation? 


a t 


1 » 


y = f (r,y)  , y(r^)  = y^  , 

the  Runge-Kutta  method  is  given  by 
It 

= E 
i=l 


yn+l-^n 


w . K . . 
i 1 


(4.1) 

(4.2) 


Here,  y =y(r  ) , the  w.  are  constants,  and 
n n 1 


i-1 


K.=h  f(r  + a.h  ,y  +.E,6..K.) 
1 n n inn  J=liJJ 


(4.3) 


with  h =r  . , -r  , and 
n n + 1 n 


a 1 =0 
^2  = 2/5 

ag  = (14-3/5)/16 

“4  “ ^ 


(4.4) 


$21  = «2 


= a 0-3 


" [“3^a3"“2^]/[2a2(l"2a2)] 


(A. 5) 


- “4“342“643 


(l-ai)ra7+a,  -l-(2a^-l)^1 

202  (013-02)  [6a2a3-4  (02+03)  + 3] 


(1-202) (I-02) (I-03) 


“3  (“3-“2)  C6O2O3  -4(a2+“3)+3. 


-2- +[1-2(02+03)^2  0203 
(2  03  -1) / [1202 (03-02) (I-02)] 


(A  .6) 


- (1-2  02)/[12o3  (o  -a  )(1-03)] 


^ ^[2(02+03)-3j/[l2(l-02)  (1-03)] 


Pcedtctor-corrector  method,  xhe  fourth-order  method  we 


have  used  (Hamming,  1959)  is  fully  described,  along  with  the  details 


concerning  doubling  of  step  size,  by  Ralston  (1960).  Highly  practi- 


tor-corrector routines,  which  we  have  employed,  will  be  found  in 


Anon. (1970),  One  should  be  warned,  however,  that  the  use  of  the 


subroutines  therein  described  is  highly  inadvisable  for  our  present 


purposes.  The  use  of  these  genei 


>urpose  subroutines  can  increase 


cal  details  concerning  the  combination  of  the  Runge-Kutta  and  predic- 


computation expense  by  a factor  of  10  to  100  over  that  of  the  optimized 


algorithm. 
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5.  OPTIMIZATION  OF  THE  ALTERMAN-JAROSCH-PEKERIS  FORMULATION 


The  key  to  optimizing  the  integration  is  to  apply  our 
knowledge  about  this  specific  problem  to  specify  all  the  depths^  r^^, 
at  which  a..(r)  are  to  be  evaluated.  The  evaluation  of  these 
elements  can  then  be  removed  from  the  innermost,  integration  loops 
of  the  program.  The  details  concerning  these  depths  are  contained 
in  the  preceding  section.  In  Figure  1,  the  optimized  scheme  for 
the  evaluation  of  a^^  (r^^) — for  the  solid  sedimentary  layers,  subsed^ 
mentary  crustal  layers,  and  mantle — ^is  indicated  in  outline  form. 

This  figure  shows  that  most  of  the  procedure  for  evaluating  a.  . (r  ) 
can  even  be  removed  from  within  theiJiand  £ loops:  within  the  c^ 
loop*  each  new  value  requires  only  6N+1  assignments,  6^+1  multiply 
cations,  and  IJ+1  additions;  within  the  ^ loop,  each  new  w value 
requires  only  314+1  assignments,  ^+1  additions,  and  2N  subtractions, 

where  N is  the  number  of  depths  at  which  a.  . (r,  ) must  be  evaluated. 

_JJ }L_ 

All  other  portions  of  the  element  determinations  are  performed 
external  to  these  loops.  In  Figure  2 the  same  information  is  given 
for  the  elements,  b^^ (r)  of  the  matrix  describing  the  integration 
through  the  liquid,  oceanic  layer.  Again,  all  elements  can  be 
evaluated  external  to  the  integration  routine. 

In  the  integration  procedures  themselves,  it  is  very 
important  to  form  matrix  products,  such  as  those  in  (2.01)  and 
(2.16),  in  an  explicit  manner.  This  permits  full  use  to  be  made 
of  the  many  zero  elements,  and  those  that  are  independent  of  £, 
or  are  equal  to  another  matrix  element.  For  example,  the  basic 
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AJP  matrix  multiplication  for  solid  layers  is  illustrated  in 
Figure  3. 


In  reports  on  our  earlier  work  with  Love-wave  dispersion 
and  dispersion-attenuation  computations,  for  both  flat  and 
spherical  structures,  it  was  possible  to  give  the  short,  key, 
FORTRAN  program  segments  (Figure  2,  Schwab  and  Knopoff,  1972; 
Figures  4 and  5,  Schwab  and  Knopoff,  1973).  These  together 
with  descriptions  of  the  root-bracketing  and  root-refining  proce- 
dures, completely  specify  the  optimization  when  the  mult i- , homoge- 
neous-layer approximation  is  employed.  When  this  approximation  is 
used  with  Rayleigh  waves  on  flat  structures,  the  optimization  can 
be  specified  in  the  same  manner  (Figures  11,  12  and  13,  Schwab  and 
Knopoff,  1972).  When  employing  the  method  of  direct  integration 
of  the  equations  of  motion,  it  is  not  possible  to  exhibit  the 
complete  program  optimizati  'n  in  this  compact,  simple  manner  for 
Rayleigh-wave  dispersion  on  a spherical,  radially  heterogeneous, 
gravitating  earth.  However,  it  possible  to  present  the  most 
important  part  of  the  algorithm  as  a relatively  compact  program 


segment.  This  is  given  in  Figure  4a, which  illustrates  the  predic- 


tor-corrector method  we  have  applied  to  the  integration  from  below 

the  Moho  to  the  selected  value  of  r ; the  automatic  doubling  of 

o^ 

integration  step  size  is  included  in  the  segment.  Most  of  the 
computation  time  is  spent  in  this  program  segment,  which  is  entered 
with : 

COEFFl  = (4/3)  H 
C0EFF2  = 3H 
C0EFF6  = (121/36)  H 
H - -25/16; 

the  indices  for  the  successive  integration  step-size  regions  are 
given  in  Table  1;  for  the  details  concerning  B(I,J),  see  the 
description  of  subroutine  DHPCL  (Anon.,  1970).  In  this  type  of 

programming  there  are  important,  machine-dependent  considerations: 


the  manner  in  which  the  6x6  matrix  elements  are  stored  in 
memory,  and  the  manner  in  which  indices  are 

handled  in  the  program  segment  given  in  Figure  4a.  In  fact, 
the  indices  ITPl,  ITP3,  ITP8,  ITP9 , ITPIO  which  are  used  for 
compactness  in  Figure  4,  actually  slow  computation  speed  on 
the  IBM  360/91;  these  subscripts  are  best  used  in  explicit 
form  IT+1 , IT+3,  IT+8,  IT+9,  IT+10.  The  key  program  segment  is 
given  in  the  form  shown  in  Figure  4a  for  two  reasons:  (1) 

to  illustrate  the  logic  as  clearly  and  simply  as  possible, 
and  (2)  to  provide  an  illustrative  example  of  the  importance 
of  handling  subscripting  and  storage  in  the  manner  most 
appropriate  for  a given  machine.  The  time  required  to 
execute  DO-loop  170  once,  specifies  the  necessary  time  to 
execute  one  integration  step  for  each  of  the  three  vectors 
; thus,  to  perform  one  integration  step  in  forming  the 
dispersion  function,  DO-loop  170  must  be  executed  three 
times.  The  time  for  one  integration  step  in  forming  ^ is 
termed  the  characteristic  time  T^,  which  we  use  to  illustrate 
the  importance  of  correct  subscripting  and  storage.  The 
characteristic  time  for  the  segment  in  Figure  4a  is 
489  psec/step/iteration . By  simply  reversing  the  order  of 
the  subscripts  of  B(I,J),  this  time  is  improved  by 
92  Msec /s t ep / i t era t ion ; if  ITPl,  etc.  are  used  explicitly 
as  IT+1,  etc.,  is  decreased  still  further  by  an  amount  of 

67  psec /s tep/ i t er at  ion ; and  if  a are  stored  more  efficiently, 

■ 

still  another  44  psec/step/iteration  can  be  saved,  bringing 
down  to  286  ps  e c / s t ep  / i t er  a t ion  . DO-loop  170,  in  a form 
incorporating  the  above  Improvements,  is  shown  in  Figure  4b. 
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It  should  be  understood  that  286  ysec/step/iteration  is  a 
lowe r bound;  the  corresponding  effective  characteristic  time 
given  in  (5.02)  --  must  reflect  time  spent  in  other  parts  of 
the  program  than  DO-loop  170,  e.g.,  time  spent  in  starting  the 
predictor-corrector  scheme  with  the  Runge-Kutta  procedure. 


Our  computations  have  been  performed  on  IBM  360  and  370 
installations:  a 360/91  in  Los  Angeles,  a 360/65  in  Bari,  and 

a 370/145  in  Santiago  and  in  Cosenza.  The  first  installation 
was  also  used  in  the  final  optimization  of  the  flat-structure 
Rayleigh-wave  work  (Schwab  and  Knopoff,  1972),  which  allows  us 
to  make  an  accurate  evaluation  of  the  relative  characteristic 
times  (Schwab  and  Knopoff,  1972)  for  computations  with  flat, 
non-gravitating  structures  and  with  spherical,  gravitating 
models.  In  the  former  case,  this  time  is 


FLAT  ^ 


RAYLEIGH 


= 110 


ysec/layer/iteration 


(5.01) 


(which  corresponds  to  Knopoff's  method  applied  to  a sequence 
of  homogeneous  layers),  and  in  the  latter  case. 


SPHERICAL 


RAYLEIGH 


=336  ysec/step/iteration  (5.02) 


(which  corresponds  to  the  optimization  of  the  AJP  formulation 
herein  described).  Again,  these  characteristic  times  were 
measured  on  the  IBM  360/91  at  UCLA.  As  we  have  noted  above. 
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an  integration  "step"  can  be  considered  nearly  equivalent  to 
a "layer"  in  computations  based  on  the  homogeneous- lay er 
approximation.  Also,  to  the  accuracy  possible  in  this  type  of 
comparison,  the  "iterations"  required  in  the  two  cases  (see 
Schwab  and  Knopoff  (1972)  for  details  concerning  iteration 
procedures)  can  be  considered  equivalent.  Thus,  the  relative 
efficiences  of  the  two  types  of  Rayleigh-wave  dispersion 
computation  can  be  evaluated  by  simple  comparison  of  their 
characteristic  times,  and  we  find  that  the  inclusion  of 
sphericity  and  gravitation  triples  actual  computation  time. 

Thus,  the  time  required  to  integrate  each  of  the  three  vectors 
over  depth  in  the  spherical,  gravitating  case,  is  the  same  as 
that  required  to  carry  out  the  analogous  operation  --  the 
formation  of  the  matrix  product  --  for  the  flat,  non-gravitating  case 
To  obtain  a valid  comparison  of  the  direct-integration 
method,  with  the  homogeneous- lay er  technique,  clearly  we  should 
not  include  gravity  in  the  former  method.  The  removal  of 
gravity  reduces  the  vectors  X^,  that  must  be  integrated,  from 
three  to  two,  and  the  number  of  elementary  operations  (multi- 
plications and  additions)  in  (2.01)  from  34  to  23.  Thus,  the 
ratio  of  computation  times  for  non-gravitating  and  gravitating 
spheres,  when  treated  with  the  direct-integration  method,  is 
approximately  (2x23)  / (3x34);  or,  for  the  non-gravitating  case. 


SPHERICAL  ^ RAYLEIGH"^^^  y se c / s t ep / i t e r a t ion 


(5.03) 


Thus  the  direct  integration  method,  for  a non-gravitating 
sphere,  is  only  36  percent  slower  than  the  optimized  com- 
putations for  a flat,  non-gravitating  structure,  where  the 
latter  is  treated  with  the  homogeneous-layer  approximation. 

From  this  result  we  are  led  to  conclude  that  attempts 
to  devise  Rayleigh-wave  transformation  techniques  --  which 
have  hitherto  been  concentrated  on  curvature  corrections  to 
permit  spherical,  gravitating  structures  to  be  treated  with 
algorithms  for  flat,  non-gravitating  models  --  might  also 
be  directed  toward  gravity  corrections  that  would  allow  one 
to  use  programs  for  spherical,  non-gravitating  structures  to 
handle  the  spherical,  gravitating  case. 

A final  Improvement  in  the  lower  bounds  of  the  AJP 
characteristic  times  is  possible.  The  limiting  bound  for 
the  gravitating  case  was  obtained  by  including  the  Integration 
of  all  three  X^,  simultaneously,  within  DO-loop  170;  the  re- 
sult was  266  y s e c / s t ep / i t e r a t ion . If  the  two  vectors  of  the 
non-gravitating  case  are  handled  simultaneously  within  the 
loop  (see  Figure  4c),  the  lower  bound  of  jr  becomes 
143  y sec/ s tep / i ter a t ion . 

It  is  obviously  of  interest,  to  those  involved  in  surface- 
wave  computations,  to  have  some  idea  of  the  relative  speeds  of 
such  computations  for  the  various  computers  currently  in  use  at 
the  larger  Installations.  In  Table  2 we  present  the  results  of 
a first  effort  to  summarize  this  information  for  the  computers 
that  are  available  to  us  for  testing.  The  test  routine  we 
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have  employed  is  that  given  by  Schwab  and  Knopoff  (Figure  2, 

1972).  A 100-layer  structure  was  used;  the  program  segment 

was  enclosed  within  a DO-loop  that  was  executed  either  100  or 

1000  times;  case  1 involved  tests  with  c<3  > and  case  2,  tests 

— m 

with  c>3  • Results  from  a more  extensive  set  of  relative  computer- 
— m 

time  tests,  for  both  Love-  and  Rayleigh-wave  dispersion  com- 
putations are  given  by  Porter  e_t  a_l.  (1977). 


6. 


EXISTENCE  OF  SOLUTIONS  AS  A FUNCTION  OF  NUMERICAL  AND 


ALGORITHMIC  PROCEDURES 

On  IBM  360  equipment,  large-scale  numerical  work  is 
routinely  carried  out  in  double  precision:  about  16  decimal 

digits.  Except  where  indicated  otherwise,  this  precision  was 
used  to  investigate  the  existence  criteria  for  solutions  from 
our  optimization  of  the  basic  AJP  formulation. 

Our  testing  procedure  followed  the  lines  of  the 
layer-reduction  experiments  described  by  Schwab  and  Knopoff 
(19  72):  At  each  of  a set  of  periods,  jc  is  computed  for  a com- 
plete range  of  terminating  values,  ro»  for  the  integration.  At 
each,  fixed  period,  by  comparing  the  values  of  c^  as  a function 
of  ij)^,  the  range  of  r_o^,  over  which  is  stable  to  4 significant 
figures  is  immediately  evident.  In  terms  of  rp  and  period,  our 
results  for  an  oceanic,  and  a continental  shield  structure 
(Figure  5)  are  given  in  Figure  6;  the  fundamental  and  first 
seven  higher  modes  are  treated  in  each  case. 

The  results  are  similar  to  those  previously  obtained  for 
Rayleigh  waves  when  the  homogeneous-layer  approximation  is 
employed  (Figures  14,  15,  17,  Schwab  and  Knopoff,  1972):  At 
each  period,  a certain  minimum  amount  of  structure  (maximum  r p ) 
must  be  retained  to  ensure  4 significant  figures  in  This 
maximum  value  of  rjj^  is  a physical  limitation.  For  the  mode  and 
period  of  interest,  there  is  significant  energy  content  down 
to  a depth  of  a- r p , and  the  structure  above  this  point  thus 


affects  all  4 figures  of  £. 
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In  the  initial  program  testing,  it  is  useful  if  one  can 
integrate  to  any  depth,  irregardless  of  whether  this  results  in  loss 
of  accuracy  in  the  computed  phase  velocity,  or  in  the  expenditure  of 
more  computer  time  than  if  the  integration  were  terminated  at  the 
minimum  acceptable  depth  for  the  desired  accuracy  in  the  phase  velo- 
city. The  use  of  a very  large  value  of  r^,  or  more  precisely,  a 
large  number  of  wavelengths  of  structure,  will  result  in  overflow; 
thus  a simple,  temporary  solution  to  this  problem  is  useful.  Such  a 
solution  is  the  simple  extension  of  the  normalization  technique 
described  by  Schwab  and  Knopoff  (1970;  1972).  The  application  of 
normalization  to  the  d ir ec t - in t e gr a t io n procedure  is  quite  simple; 
it  is  not  necessary  to  begin  normalization  until  the  application  of  the 

predictor-corrector  method  has  begun  in  the  mantle.  To  normalize, 
all  one  need  do  is  determine  the  maximum  of  the  absolute  values  of 
y^  at  the  end  of  each  integration  step;  one  then  divides  all 
y^(r^)  and  y ^ ( r ^ ) by  this  value,  where  ^ are  the  seven  positions 

at  which  y^  and  y.  must  be  specified  so  as  to  permit  the  next  step 

of  the  fourth-order  predictor-corrector  method.  Seven,  rather  than 
four  r^  are  required  to  allow  the  automatic  doubling  of  step  size 
when  certain  depths  are  reached. 

For  ease  of  reference,  a normalization  scheme 
which  is  appropriate  for  the  program  segment  in  Figure  4a,  is 
given  in  Figure  7.  Two  warnings:  (1)  If  only  sparing  use  of 

normalization  is  planned,  or  it  is  to  be  invoked  from  an  IF 
statement,  the  segment  in  Figure  7 will  be  satisfactory;  how- 
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ever,  if  large-scale  use  is  envisioned,  efficiency  requires  the 
inclusion  of  normalization  directly  within  the  coding  in  Figure  4a. 
(2)  Unless  absolutely  necessary,  normalization  should  not  be  in- 
cluded in  these  computations;  it  can  result  in  a very  significant 
increase  in  computation  time. 

Our  numerical  tests  of  the  overflow  problem  were  performed 

with  the  average  (oceanic)  earth  structure  given  by  Wiggins  (1968). 

The  results  are  given  in  Table  3.  For  IBM  360  equipment,  when  using 

7 0 8 0 

double-precision  computation,  overflow  occurs  when  1^1  = 10  to  10 

Returning  to  Figure  6,  for  routine  use  of  the  AJP 
formulation,  it  is  important  to  integrate  only  down  to  slightly 
below  the  position  the  Tj^  lines  in  the  figure,  and  to  thereby 
minimize  computation  time  and  expense.  For  this  purpose  we  have 
devised  empirical  "laws"  for  determining  the  maximum  depth 

H = a - r 0 

to  which  the  integration  must  be  performed  if  the  resulting 
phase  velocity  is  to  be  valid  to  4 significant  figures.  The 
data  for  these  determinations  are  collected  in  Figure  8.  The 
"laws"  specifying  the  number  of  wavelengths  of  structure  to  be 
retained,  if  4 significant  figures  are  desired  in  the  computed 
phase  velocities,  are 


H/A 

= 7-c 

for 

mode  0 

(6.01) 

H/A 

-9.5-7 

for 

mode  1 

(6.02) 

H/A 

- [11  + ^ (M-2)]  - c 

10 

for 

modes  2-7 

(6.03) 

where  \ is  the  wavelength,  and  M is  the  mode  number. 
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From  the  results  shown  In  Figures  6 and  8 we  also 
have  the  maximum  periods  and  values  of  c which  can  be  determined 
without  taking  the  structure  of  the  core  into  consideration 
(other  than  to  determine  the  values  of  g (r) ) . These  results 
are  shown  in  Figure  9,  along  with  the  corresponding  values  of  the 
minimum  order  number  £. 

Relative  to  the  computation  of  theoretical  seis- 
mograms, the  combination  of  the  results  in  Figure  9 and  those 
given  by  Schwab  and  Kauael  (Section  5,  1976),  indicate  a potentially 

useful  conclusion:  (1)  When  li>SL  . , only  the  crust-mantle 

— min 

system  need  be  used  in  the  computations;  c can  be  computed  at 
specified,  equally-spaced  frequencies  and  inverse  Fourier  trans- 
formation can  be  used  to  calculate  the  theoretical  seismogram  for 
this  range  of  periods;  and  the  first  term  of  the  asymptotic  ex- 
pansions for  and  can  be  used  (possibly  corrected  by 

automatic  numerical  interpolation  from  the  data  in  Figures  2 and 
3 of  Kausel  and  Schwab  (1976)).  (2)  When  ^<IL  the  core 

must  be  included  in  the  computations;  w should  be  computed  at 
integral  values  of  and  summation  over  ^ should  replace 
Fourier  synthesis;  and  the  exact,  integral-^  expressions  should 
be  used  to  evaluate  the  associated  Legendre  functions. 


A point  of  considerable  interest  to  those  involved 
in  the  actual  computation  of  surface  wave  dispersion,  is  whether 
or  not  a difficulty  analogous  to  the  Thomson- Haskell  "loss-of- 
precision"  problem(Schwab  and  Knopoff,  1970;  Schwab,  1970)  is 
intrinsic  to  the  AJP  formulation.  In  computations  based  on 
the  homogeneous-layer  approximation,  when  the  original  version 
(Haskell,  1953)  of  the  Thomson-Haskell  formulation  for  Rayleigh 
waves  is  used,  this  problem  can  cause  serious  difficulties  if 
the  computer  is  employed  in  a low-precision  mode.  To  test  for 
an  analog  to  this  "loss-of-preclsion"  problem,  in  our  optimization 
of  the  AJP  formulation,  we  simulated  single-precision  (about  6 
decimal  digits)  computation  by  replacing  DO-loop  160  (Figure  4a) 
in  our  double-precision  program,  with  the  program  segment  shown 
in  Figure  10.  The  function  SNGL  accepts  a double-precision 
argument,  andreturns  the  single-precision  equivalent. 

The  results  of  our  single-precision  tests  are 
similar  to  those  from  the  original  Thomson-Haskell  formulation 
(Figure  2,  Schwab  and  Knopoff,  1970),  and  are  illustrated  in 
Figure  11.  In  the  range  shown,  our  results  indicate  that 

there  is  no  problem  with  modes  0-4  in  double-precision  computations; 
but  when  computations  are  reduced  to  single  precision,  the  loss-of- 
preclslon  problem  is  clearly  illustrated.  In  the  latter  case, 
there  is  seen  to  exist  a minimum  value  of  rj^,  below  which  we 
cannot  go  and  still  retain  a given  accuracy  in  the  computed  phase 
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velocities.  Thus  the  AJP  formulation  does  Indeed  exhibit  the 
analog  to  the  Thomson-Haskell  loss-o£-pr eels  Ion  problem. 

Since  practical  work  with  dispersion  computations 
on  IBM  360  equipment  Is  routinely  carried  out  In  double  precision, 

It  Is  important  to  estimate  the  numerical  limitations  imposed  by 
the  loss-of -precis ion  problem  in  this  computational  mode.  The 
results  of  our  tests  at  50  and  25  seconds  are  shown  In  Figures 
12  and  13.  Less  extensive  tests  were  also  carried  out  at  a 
period  of  65  seconds.  At  a given  period,  the  right-most  point  of 
each  of  the  smoothed  curves  was  used  to  determine  the  maximum 
accuracy  possible  for  each  mode.  This  Information  was  then 
collected  in  the  Figures  14  and  15.  Although  the  data  is 
necessarily  sparse,  due  to  the  expense  of  this  type  of  experiment, 
the  results  are  strikingly  clear:  For  a fixed  period,  as  we  go 
to  higher  and  higher  mode  numbers,  the  attainable  accuracy  In 
becomes  less  and  less;  for  a fixed  accuracy  In  c_,  as  we  go 
to  shorter  and  shorter  periods,  the  maximum  mode  number  that 
can  be  successfully  treated  becomes  smaller  and  smaller. 

In  the  near  future  we  Intend  to  present  the  details  of 
our  numerical  analysis  of  the  various  possible  methods  for  dealing 
with  the  loss-of-precls ion  problem.  See,  for  example,  Wiggins  (1968), 
Nolet  (Appendix  A,  1976),  Nelgauz  and  Skadlnskaya  (1972),  Gilbert 
and  Backus  (1966),  and  Takeuchl  and  Salto  (Section  II. D. 4,  1972). 


£ 
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7.  TERMINATING  BOUNDARY  CONDITIONS 


In  our  present  programming  efforts  we  are  concentrating 
on  the  computation  of  phase  velocities  at  arbitrarily  specified 
(equally-spaced)  frequencies.  To  keep  our  algorithms  as  simple 
as  possible,  and  to  avoid  the  difficulties  involved  in  the 
evaluation  of  spherical  Bessel  functions  of  non-integral  order, 
we  have  chosen  to  terminate  our  integrations  at  depth  with  free, 
or  rigid,  terminating  boundary  conditions.  For  completeness,  we 
also  include  the  technique  for:  (1)  terminating  the  integration 

within  the  mantle  by  applying  terminating  boundary  conditions  for 
a gravitating,  homogeneous,  solid  sphere  below  r^  and,  (2)  termi^ 
nating  at  the  mantle-core  boundary  by  applying  the  conditions  for  a homo- 
geneous liquid  sphere  below  r^^ 


In  the  former  case,  just  above  r^  we  have 


) = yi(a)Xi(ro)  + y3(b)X2(ro)  + y5(a)X3(ro)  (7.01) 


where 


a for  a continental  structure 


(7.02) 


r^  for  an  oceanic  structure 
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For  the  homogeneous,  gravitating,  solid  sphere  below  r^^  , there 
are  three  classes  of  solutions:  Yi(r),  YsCr);  thus  just 
be  low  , 

Y_(ro)  = D Yi(ro)  + E Y2(ro)  + M Yalro),  (7-03) 
where  D,  E,  and  M are  undetermined  constants.  Applying  the  boundary 


conditions 

of  continuity  of  y^ 

at  rp , we 

obtain 

[Xi (ro)]  1 

[_X2(ro)]  1 

[X3(ro)J 1 

-[Yi(ro)]i 

"LX2(>^o)J  1 

-[Y3(’^o)] 

[Xi(ro)]2 

[X2(ro)]2 

[X3(*^o)]2 

CM 

O 

U 

-&2(ro)32 

-[Y3(ro)] 

[Xi (ro)  J 3 

[X2(ro)l 3 

[Xa  (^^0)3  3 

-[Yi(ro)j3 

-LY2(’^o)J  3 

-[Y3(ro)] 

[Xi(ro)]4 

[X2(ro)]4 

[X3(ro)]4 

-[^1  (ro)j4 

-[Y2(ro)]4 

-[^3(^0)] 

[Xi (ro)  1 5 

[X2(ro)]5 

[X3(>^0)]5 

-[Yl(ro)]5 

-[^2(^0)]  5 

-[Y3(ro)] 

1 [Xi(ro)l6 

[Xa (ro)]6 

[X3(ro)]6 

-[Yi(ro)]6 

-[Y2(J^o)]6 

-[Y3(ro)] 

(7.04) 


NW  = 0 

and  the  dispersion  function  has  the  form 


(7.05) 


(7.06) 


F(aj,c)  = det(N)  . 
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The  components  of  Y^(r)  are  given  in  convenient  form  by  Takeuchi 

and  Saito  (1972):  YjCr)  bytheir  equations  (98)^with  the  negative 

sign  in  (99);  Y2(r)  by  (98),  with  the  positive  sign  in  (99);  and  Y^(r) 

by  their  equations  (100).  Note  that  their  definition  of  y.  differs 

6 

slightl/  from  that  used  here. 

When  the  structure  used  to  form  the  dispersion  function  is 
terminated  at  the  mantle-core  boundary  by  the  conditions  for  a 
homogeneous  liquid  sphere  below  r^  , above  r^^  we  have  (7.01);  below, 

Y_(ro)  = P Yi(ro)  + Q Y3(ro)  , (7.07) 

where  ^ and  ^ are  undetermined  constants, 
and  Y^  in  (7.07)  have  the  form 

Y.(ro) 

Again,  Takeuchi  and  Saito  (1972)  give  the  form  of  these  vector 
components  for  the  gravitating,  homogeneous,  liquid  sphere.  From 
the  conditions  of  continuity  of  yj,  y2 . ys.  ye 
vanishing  of  y4(ro),  we  have 
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~ n 

— 

[?^l(*^o)-ll  [X2(>^o)ll  [X3(j^o)]l  “[Yi(ro)]l  "&3(ro)]i 

yi  (a) 

0 

[,Xi(ro)j2  Lx2(ro)]2  CX3(ro)]2  -[Yi(ro)]2  -&3(>^o)]2 

y3  (b) 

0 

[Xi(ro)_]‘4  [X2(ro)]4  [X3(^o)]4  0 0 

y5  (a) 

0 

lXl(’^o)J5  [X2(^^o)]5  [X3(>^0)]5  “&l(l'o)l5  -&3(’^o)35 

P 

0 

[Xi(ro)]6  [^2^’^0^36  p3^’^0^]6  ~p3^*^0^J5 

Q 

0 

or 


R S = 0 , 


(7.09) 


(7.10) 


and  the  dispersion  function  takes  the  form 


F(a),c)  = det(R)  . (7.11) 

To  obtain  the  group  velocity  we  still  employ  (3.03)  and  (3.04). 

The  forms  of  and  F^ which  result  from  (7.06),  or  (7.11),  can  be  obtained  by 
analogy  with  the  way  in  which  (3.05)  is  obtained  from  (2.10).  In  the  present 
case,  however,  the  analog  of  (3.05)  will  comprise  the  sum  of  six,  sixth-order 
determinants  when  (7.06)  is  used  to  form  the  dispersion  function,  and  the  sum 
of  five,  fifth-order  determinants  when  (7.11)  is  used. 

As  would  be  expected,  our  numerical  tests  of  a 
terminating  solid  sphere  show  that,  to  obtain  a given  accuracy 
in  c , less  structure  must  be  retained  in  this  case  than  when 
using  terminating  rigid  or  free  boundaries  within  the  solid  mantle. 
A complete  set  of  tests,  comparable  to  those  in  Figure  12,  was 
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performed  for  a period  of  50  seconds;  the  results  showed  that 
the  increase  in  the  maximum  rj^  values  was  surprisingly  independ- 

l ent  of  mode  number  and  o : (300+25)  km,  or  in  most  cases, 

I -n  ~ 

about  24  fewer  integration  steps  when  terminating  with  a solid 
sphere.  Thus  the  use  of  the  "correct"  terminating  boundary 
condition  appears  to  be  important  only  for  the  lowest  (radial) 

I mode  numbers.  Our  tests  with  mode  7 at  25  seconds,  show  the 

i 

increase  in  r p to  be  about  (110+10)  km.  Thus,  from  our  limited 
number  of  tests,  it  appears  that  Sto/tq  Is  roughtly  constant  with 
about  the  value  0.14+0.02  ; where  rj^  is  maximum  value  required 
by  the  structural  limitation  when  the  condition  of  a rigid 

I 

terminating  boundary  is  employed,  and  6r p is  the  increase  possible 
in  when  one  then  employs  the  condition  of  a terminating  solid 
sphere . 

Tests  of  the  loss-of-precision  problem  were  also  per- 
formed with  the  "correct"  boundary  condition  at  depth.  Again,  a com- 
plete set  of  tests  was  carried  out  at  50  seconds,  and  mode  7 was 
tested  at  a period  of  25  seconds.  The  right-hand  extremes  of  the 
analogs  of  the  smooth  curves  in  Figures  12  and  13,  occurred  at  the 

I same  depths  in  these  new  tests;  thus,  as  a result  of  the  Increased 

r 0 values  of  the  upper  portions  of  the  curves  when  solid-sphere 
termination  is  used,  the  maximum  accuracy  for  any  given  mode  is 
significantly  improved  by  using  this  type  of  boundary  at  depth. 
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8. 


CONCLUSIONS 


An  analysis  of  direct,  Rayleigh-wave  dispersion  com- 
putations on  a spherical,  gravitating  earth  has  been  performed 
using  the  Alterman-Jarosch-Pekeris  (1959)  formulation. 

1.  No  difficulty  was  encountered  when  we  reversed  the 
usual  procedure,  for  practical  purposes,  and  computed  phase 
velocities  (or  polar  order  numbers)  at  specified  periods. 

2.  Integration  from  the  free  surface  downward,  again 
reversing  the  "standard"  procedure,  resulted  in  no  unexpected 
difficulties.  In  fact,  this  procedure  much  simplified  the 
specification  of  the  algorithm  for  integrating  the  system  of 
differential  equations  to  obtain  phase-velocity  dispersion. 

This  procedure  makes  the  generalization  from  the  algorithm  for 
continental,  to  oceanic  structures  relatively  trivial;  also, 
it  makes  it  possible  to  develop  direct  algorithms  for  obtaining 
group  velocities  for  the  two  types  of  structures.  The  usual 
variational  techniques  are  not  required  for  this  latter  purpose. 

3.  Our  optimization  of  the  AJP  formulation  is  based  on 

removing  all  function  evaluations  from  the  innermost,  Integration 

(over  r^)  loops  of  the  program.  In  fact,  most  of  the  evaluation 

procedure  for  a (r  ) can  even  be  removed  from  the  phase  velocity 
. j-J  i 

and  frequency  loops.  This  optimization  of  the  AJP  formulation 
yields  a characteristic  time  of  336  yseconds /integration  step/ 
Iteration  for  spherical,  gravitating  structures,  and  a time  of 
151  yseconds/lntegratlon  step/iteration  for  spherical. 
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non-gravitating  structures.  These  characteristic  times,  for 
the  AJP  direct-integration  procedure,  provide  a basis  of  com- 
parison with  the  homogeneous-layer  approximation  for  a non- 
gravitating flat  structure,  which  has  a characteristic  time 
of  110  y se c / lay e r / i t e r a t ion  (Knopoff's  method,  Schwab  and 
Knopoff,  1972).  The  lower  bounds  of  the  characteristic 
times,  for  our  final  optimizations  of  the  AJP  formulation,  are 
266  y sec /s t ep/ i terat ion  for  the  gravitating  case  and 
143  ysec/step/iteration  when  gravity  is  not  included.  All  of 
the  above  times  apply  to  the  IBM  360/91  computer  at  UCLA. 

4.  Our  results  here ^ combined  with  those  of  Schwab 

and  Knopoff  (1972),  indicate  that  an  integration  "step"  (in 

the  AJP  procedure)  can  be  considered  nearly  equivalent  to  a "layer" 

in  computations  based  on  the  homogeneous- lay er  approximation. 

Also,  to  the  accuracy  possible  in  this  type  of  comparison,  the 
"iterations"  required  in  the  two  techniques  (see  Schwab  and 
Knopoff  (1972)  for  details)  can  be  considered  equivalent.  Thus 
the  relative  efficiencies  of  the  two  types  of  Rayleigh-wave  dis- 
persion computations  can  be  evaluated  by  simple  comparison  of  the 
above  characteristic  times.  The  fact  that  approximately  the  same 
number  and  sizes  of  "steps"  must  be  used  in  the  direct-integration 
procedure,  as  number  and  sizes  of  "layers"  in  the  homogeneous- 
layer  approximation,  means  that  the  usual  assumption  that  the 
former  method  does  a better  job  of  treating  continuous  parameter- 
depth  distributions,  appears  to  be  invalid. 

5.  The  overflow  problem  in  the  AJP  formulation  can  be 


controlled  by  simple  normalization.  Program  segments  are  given 
which  describe  the  procedure  explicitly. 
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6.  A loss-of-precision  problem  appears  to  be  Intrinsic 
to  the  AJP  formulation.  Results  of  this  problem:  For  a fixed 
period,  as  mode  number  increases,  the  attainable  accuracy  in 
the  phase  velocity  decreases;  for  a fixed  accuracy  in  the  phase 
velocity,  as  period  decreases  the  maximum  mode  number  that  can  be 
treated  successfully  decreases. 
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Table  1.  Constants  for  integration  through  successive  step-size  regions 
illustrated  in  program  segment  given  in  Figure  4.  Constants 


correspond 

to  4 significant  figures  in  computed  phase  velocity 
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N1(I) 

N2(I) 

Integration 
step  size 
(km) 
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11 

-1.5625 
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12 

16 

-3.1250 

3 

17 

21 

-6.2500 
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22 

* 

-12.5000 

* N2(4) 

is  specified 

so  as 

to  allow  integration  to  proceed  to  the 

deepest 

point  within 

the  solid  mantle,  while  maintaining  a step 

size  of 

-12.5  km.  NEND  is 

determined,  at  each  period,  by  the  input 

value  of  r,;  it  must  satisfy  NEND  ^ N2(4). 
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Table  3.  Results  of  numerical  tests  of  the  overflow 


problem  when  normalization  is  not  Included  in  our 


optimization  of  the  basic  AJP  formulation.  An 


average  (oceanic)  earth  structure  (Wiggins,  1968), 

and  a period 

of  50  seconds,  were  used  in 

the  tests. 

The  value  of 

ro  is  the  maximum  at  which 

overflow 

occurs;  H/X 

is  the  corresponding  number 

of  wave- 

lengths  of  structure,  from  the  surface  of  the  earth 

down  to  r = r 0 . 

These,  and  larger  values 

of  H/X,  yield 

overflow , 
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FIGURE  CAPTIONS 


Figure  1.  Schematic  representation  of  optimized  scheme  for 

evaluating  the  matrix  elements  a (r.)  In  the 

^ 

treatment  of  the  solid  sedimentary  layers,  the 
subsedimentary  crustal  layers,  and  the  mantle. 

The  quantities  \ and  ^ are  Lame's  constants,  £ 

Is  the  gravitational  constant,  and  N Is  the  number 
of  depths  at  which  a^*^  must  be  evaluated. 

Figure  2.  Schematic  representation  of  optimized  scheme  for 

evaluating  the  matrix  elements  b (r.)  in  the 

!S_ 

treatment  of  the  homogeneous  oceanic  (liquid)  layer. 

The  quantity  a is  the  compressional-wave  velocity, 

£ is  the  density,  g (r)  is  acceleration  due  to  gravity, 

and  M Is  the  number  of  depths  at  which  b (r.)  must 

^3  !S_ 

be  evaluated. 

Figure  3.  (a)  FORTRAN  IV  program  segment  for  the  basic  matrix 

multiplication  for  our  optimization  of  the  AJP 
formulation  for  solid  layers;  (b)  symbolic  repre- 
sentation of  (a),  which  Is  used  In  Figure  4;  and 
(c)  definition  of  one-dlmenslonal  array  used  In  (a) 
to  represent  6x6  matrix  In  (2.01).  The  Integer  IPT 
>■  Is  the  Index  specifying  the  value  of  £,  and  A1 

through  A36  are  dimensioned  to  300. 

Figure  4a.  FORTRAN  IV  program  segment  In  which  the  predictor- 

corrector  portion  of  the  Integration  from  below  the  Moho 
to  rg  is  handled;  most  of  the  computation  time  Is  spent 
In  this  segment.  The  boxed  segments  refer  to  the  basic 


Figure  4b.  FORTRAN  IV  program  segment  demonstrating  sub- 
scripting and  storage  improvements,  relative 
to  the  segment  in  Figure  4a,  that  are  required 
to  optimize  computation  time  on  an  IBM  360 
computer . 


Figure  4c.  FORTRAN  IV  program  segment  illustrating  our 

final  optimization  for  the  non-gravitating  case 


Figure  5. 


Figure  6. 


Figure  7 . 


Figure  g. 


Figure  9. 


Oceanic  and  continental  (shield)  models  used  in 
program  testing. 

Values  of  the  depth  at  which  integration  is 

terminated,  which  yield  4-significant-figure 
accuracy  in  the  computed  values  of  £ with  the 
optimized  version  of  the  basic  AJP  formulation. 

At  each  period,  4-figure  accuracy  is  attained 
only  if  to  is  specified  to  be  smaller  than  the 
indicated  curve  . 

Normalization  scheme  appropriate  to  program  segment 
in  Figure  4a.  The  procedure  should  be  included 
between  statement  numbers  160  and  170  in  Figure  4a. 
See  text  for  warnings  concerning  loss  of  efficiency 
when  normalization  is  employed. 


Results  of  tests  for  determining  general,  multimode 
"laws"  for  specifying  the  required  values  of  r^  to 
use  when  computing  phase-velocity  dispersion  for 
mantle  Rayleigh  waves,  if  c is  desired  to  4-figure 
accuracy. 


Maximum  periods,  maximum  phase  velocities,  and 
minimum  order  numbers,  that  can  be  used  if 

c is  desired  to  4-figure  accuracy,  when  the  inte- 
gration is  limited  to  the  mantle  and  the  core  is 
excluded  from  the  computations  (other  than  for  use 
in  the  determination  of  g (r ) ) . 


Figure  10.  FORTRAN  IV  Program  segment  used  to  simulate 
single-precision  computations  when  using  our 
double-precision  optimization  on  the  basic  AJP 
formulation.  This  program  segment  is  used  to 
replace  DO-loop  160  in  Figure  4a. 

Figure  11.  Test  of  the  effect  of  reducing  IBM  360  com- 
putations from  double  (about  16  decimal  digits) 
to  single  precision  (about  6 decimal  digits), 
while  keeping  period  fixed  at  50  seconds.  For 
a given  mode,  in  order  to  obtain  a_  significant 
figures  in  the  computed  phase  velocity,  r_o^  must 
not  exceed  the  value  given  by  the  upper  portion 
of  the  dashed  line  (structural  limitation)  nor 
fall  below  the  lower  portion  (loss-of-precision 
limitation) . 

Figure  12.  IBM  360  double-precision  tests  of  loss-of-precision 
problem  at  a period  of  50  seconds.  At  left,  raw 
results  of  structure  reduction  experiments;  at  right, 
smoothed  curves  for  each  mode.  Latter  curves  are 
drawn  such  that  all  data  points,  for  a given  mode, 
fall  to  the  right  of  corresponding  curve. 


Figure  13.  IBM  360  double-precision  tests  of  loss-of-precision 
problem  at  a period  of  25  seconds. 


Figure  14.  Maximum  possible  mode  number,  n , for  which 

max 

£ significant  figures  can  be  obtained  with  our 
optimization  of  the  AJP  formulation.  This 
limitation  is  due  to  the  loss-of-precision 
p r ob lem . 


Figure  15.  Relationship  between  maximum  possible  mode  number, 

n , and  minimum  period,  for  several  values  of  a. 
max  — 

This  limitation  of  our  optimization  of  the  AJP 
formulation  is  due  to  the  loss-of-precision  problem. 


Elements  to  be  evaluated  external  to  both  o)  and  c loops: 


% 

3l  1 

a 1 2 

333 

345 

34  1 

32  2 

334 

32  6 

r 

35  1 

34  2 

34  4 

3e  6 

Auxiliary 

quantities 

to  be  evaluated 

external  to  both  o)  and 

c loops : 

2 

d2 i = 4y (3X+2p) aj  2/r 

2 

di,  3=-2y/r 

di 3=Xai 2/r 
de  3=-4ttGp 
de  5=l/r 

2 

04  3 = 4y (X+y) ai 2/r 

— ui  loop 

2 

OMEGSQ=u) 

DO  10  1=1, N 
TEMP=-a2  6 (I) ^OMEGSQ 
321 (I)=d2i (I)-TEMP 
10  £43 (I)=d43 (I) -TEMP 
r-  C loop 

0RDER=£  («.+!) 

DO  20  1=1, N 

ai  3 (I)=di 3 (I) BORDER 

323 (1)— 341 (I) BORDER 

I 

3g 3 (I)=d6 3 (I) BORDER  I 

I 

324  (I)— 333  (I)  ORDER  ] 

365(I)=d65(I)^ORDER  f 


20  34  3 (I)=fi.3  (I)+e4  3 (I)><ORDER 


Evaluate  bsG  external  to  both  u)  and 
quantities  to  be  evaluated  external 
2 

XLAINV=l/a  p P65=P62P 


c loops 
to  both 

P2  1 = 


I lii 


2 

Pi 2=-l/pr 


2 

Ps  2=4TiG/r 
2 


Pi 5=-l/r 


loop 

2 

RHMOSQ=-pw 

2 

OMSQl."'=l/u) 
DO  10  1=1, M 


P6l=-P65g(r)  S21= 

2 

P22=-g(r)/r 
P2  5=P2  2 P 


ha  1 (I)=P2  1 (I)  ><0MSQ1N 
q2 1 (I) =S2 1 (I) +RHMOSQ 


he  1 (I)=P6 1 (I) xOMSQlN 
hi  2 (I) =Pi 2 (I) xOMSQlN 
ha  2 (I) =P2  2 (I) xOMSQlN 


h62 (I)=P62(I)^0MSQIN 
hi  5 (I)=Pi 5 (I) xOMSQIN 
has  {I)=P25 (I)xOMSQIN 
10  he  5 {I)=P6  5 (I) XOMSQIN 

r—  C loop 

ORDER=£  («,+!) 

DO  20  1=1, M 

b22  (I)=ha2 (I) xORDER 

bi  1 (I)=b6G  (I)-ba2 (I) 

ba 1 (I)=ha 1 (I) xORDER+qa i (I) 

bG  1 (I)=hG 1 (I) XORDER 

b 1 a ( I ) =h 1 a ( I ) xORDER+XLAINV 

bG2  (I)=h6  2(I)>‘ORDER 

bi  5 (I)=hi 5 (I) XORDER 

ba 5 (I) =h2 5 (I) XORDER 


Auxiliary 
u)  and  c loops: 

P2  sg (r) 

4Tig  (r)  /r 


yBAR(l)=Al(IPT)  *Y(1)+A7  (IPT)  *Y(2)+A13  (IPT)  *Y(3) 
YBAR(2)=A2  (IPT)  *Y{1)+A8  (IPT)  *Y(2)+A14  (IPT)  *Y(3)  + 

1 A20 (IPT) *Y(4) +A32 (IPT) *Y(6) 

(a)  YBAR(3)=A15(IPT) *(  Y (3) -Y (1) ) +A21 (IPT) *Y (4) 

YBAR(4)=A4 (IPT) *Y (1) +A10 (IPT) * Y(2) +A16 (IPT) *Y(3) + 
1 A22(IPT) *Y(4)+A28(IPT) *Y(5) 

YBAR(5) =A5 (IPT) *Y (1) +Y (6) 

YBAR(6)=A18(IPT) *Y ( 3) +A30 (IPT)  *Y (5) +A36 (IPT) *Y(6) 


(b) 

[ 

YBAR, 

IPT 

A1 

A7 

A13 

0 

0 

0 

A2 

A8 

A14 

A20 

0 

A32 

-A15 

0 

A15 

A21 

0 

0 

A4 

A5 

AlO 

0 

A16 

0 

A22 

0 

A28 

0 

0 

1 

0 

0 

A18 

0 

A30 

A36 

r 


C BEGIN  APPLICATION  OF  PREDICTOR-CORRECTOR  METHOD. 

DO  110  1=1,6 
110  PMNUSC(I)=0,0D+00 

HH=0 . 5D+00*1. 5625D+00 

C LOOP  OVER  REGIONS  WITH  DIFFERENT  STEP  SIZES. 

DO  180  IREG=1,NUMREG 

HH=HH+HH 

NSTART=N1 (IREG) 

NSTOP=N2 (IREG) 

NTEMP=5 

IT=0 

C LOOP  OVER  DEPTH  IN  CURRENT  STEP-SIZE  REGION. 

DO  170  N=NSTART,NSTOP 
IF(IT.EQ.4)  GO  TO  115 
IT=NTEMP-4 
ITP1=IT+1 
ITP3=IT+3 
ITP8=IT+8 
ITP9=IT+9 
ITP10=IT+10 
115  DO  120  1=1,6 
C SET  PREDICTOR  P{I). 

P(I)=B(IT,I)+COEFFl* (2.0D+00* (B (ITPIO , I) +B (ITP8 , 1)  ) 
-B (ITP9 ,1) ) 

C SET  MODIFIED  PREDICTOR  XM{I). 

120  XM{I)=P(I) -.9256198347107438D+00*PMNUSC(I) 

IPT=IPT+1 

XMBAR, IPT 

DO  130  1=1,6 
C SET  CORRECTOR  C(I) 

C{I)=.125D+00* (9.0D+00*B(ITP3,I) -B(ITP1,I) 

+COEFF2  * ( XMBAR ( I ) +2 . OD+0  0 *B ( ITPl 0 , 1 ) 
-B(ITP9,I))) 

PMNUSC(I)=P{I)-C(I) 

C SET  SOLUTION  VECTOR  AT  NTH  DEPTH. 

130  Y(I)=C(I) +.07438016528925620D+00 

*PMNUSC(I) 

IF (N.EQ.NEND)  RETURN 

C SET  DERIVATIVE  OF  SOLUTION  VECTOR  AT  NTH  DEPTH. 

YBAR,  JPT 

IF (NTEMP.GT. 7)  GO  TO  150 
NTMPP7=NTEMP+7 
DO  140  1=1,6 
B(NTEMP,I)=Y{I) 

140  B(NTMPP7,I)=YB^R(I) 

NTEMP=NTEMP+1 
GO  TO  170 
150  DO  160  1=1,6 
B(1,I)=B(2,I) 

B(2,I)=B{3,I) 

B(3,I)=B(4,I) 

■n  r /I  ry  ~t>  f t > 


I 


C 

C 


C 

C 


B(5,I)=B(6,I) 

B{6,I)=B{7,I) 

B(7,I)=Y(I) 

B(8,I)=B(9,I) 

B(9,I)=B(10,I) 

B(10,I)=B{11,I) 

B(11,I)=B(12,I) 

B(12,I)=B(13,I) 

B(13,I)=B(14,I) 

160  B(14,I)=YBAR(I) 

170  CONTINUE 

RESET  STORED  VALUES  OF  COEFFICIENTS  IN  PREPARATION  FOR 
DOUBLED  STEP  SIZE. 

COEFFl=COEFFl+COEFFl 

COEFF2=COEFF2+COEFF2 

COEFF6=COEFF6+COEFF6 

RESET  STORED  VALUES  OF  Y ( I ) , YBAR ( I ) , and  PMNUSC(I)  IN 
PREPARATION  FOR  DOUBLED  STEP  SIZE. 

DO  180  1=1,6 

PMNUSC{I)=8.962962962962963D+00* (Y(I) 

-B(1,I) ) -COEFF6* (YBAR(I) +B(8,I) 
+3.0D+00* (B(12,I) +B(10,I) ) ) 

B{2,I)=B(3,I) 

B(3,I)=B(5,I) 

B(4,I)=B(7,I) 

B(9,I)=B(10,I) 

B(10,I)=B(12,I) 

180  B{11,I)=B(14,I) 


DIMl-NSION  B(6,1A)  ,Y(6)  ,YBAR(6)  ,XM(6)  ,XMBAR(6)  ,P(6)  ,C(6)  ,PMNUSC(6)  , 
1 A(2090) 

EQUIVALENCE  (Y(1),Y1) 

EQUIVALENCE  (XM(I),XM1) 

EQUIVALENCE  (XMBARU)  .XMBARI) 

IPT=-18 


C LOOP  OVER  DEPTH  IN  CURRENT  STEP-SIZE  REGION. 

DO  170  N=NSTART,NSTOP 
IF(IT.EQ.A)  GO  TO  115 
lT=NTEMP-4 
115  DO  120  1=1,6 
C SET  PREDICTOR  P(I)  . 

P (I) =B (I , IT)+C0EFF1* (2 . OD+OO* (B (I , IT+10)+B (I , IT+8) ) -B (I , IT+9) ) 

C SET  MODIFIED  PREDICTOR  XM(I) . 

120  XM(I)=P(I)-.92561983A7107A38D+00*PMNUSC(I) 

IPT=IPT+19 

XMBAR(1)=A(IPT)*XM(1)+A(IPT+1)*XM(2)+A(IPT+2)*XM(3) 

XMBAR(2)=A(IPT+3)*XM(I)+A(IPT+4)*XM(2)+A(IPT+5)*XM(3) 

1 +A(1PT+6)*XM(4)+A(1PT+7)*XM(6) 

XMBAR(3)=A(IPT+8)*(XM(3)-XM(1))+A(IPT+9)*XM(4) 
XMBAR(4)=A(IPT+10)*X>I(1)+A(IPT+11)*XM(2)+A(IPT+12)*XM(3) 

1 +A(IPT+13)*XM(4)+A(IPT+14)*XM(5) 

XMBAR(5)=A(IPT+15)*XM(1)+XM(6) 

XMBAR(6)=A(IPT+16)*XM(3)+A(IPT+17)*XM(5)+A(IPT+18)*XM(6) 

DO  130  1=1,6 
C SET  CORRECTOR  C(I)  . 

C ( I)  = . 125D+00* (9 . 0D+00*B (I , IT+3) -B (I , IT+l)+COEFF2* (XMBAR (I)+2 . OIH-0 
1 *B(l,IT+10)-B(I,IT+9))) 

PMNUSC(I)=P(I)-C(I) 

C SET  SOLUTION  VECTOR  AT  NTH  DEPTH. 

130  Y(I)=C(I)+.07438016528925620D+00*PMNUSC(I) 

IF(N.EQ.NEND)  RETURN 

C SET  DERIVATIVE  OF  SOLUTION  VECTOR  AT  NTH  DEPTH. 

YBAR(1)=A(IPT)*Y(1)+A(1PT+1)*Y(2)+A(IPT+2)*Y(3) 

YBAR(2)=A(IPT+3)*Y(1)+A(1PT+4)*Y(2)+A(IPT+5)*Y(3) 

1 +A(1PT+6)*Y(4)+A(IPT+7)*Y(6) 

YBAR(3)=A(IPT+8)*(Y(3)-Y(1))+A(IPT+9)*Y(4) 
YBAR(4)=A(1PT+10)*Y(1)+A(IPT+11)*Y(2)+A(1PT+12)*Y(3) 

1 +A(IPT+13)*Y(4)+A(IPT+14)*Y(5) 

YBAR(5)=A(IPT+15)*Y(1)+Y(6) 

YBAR(6)=A(IPT+16)*Y(3)+A(1PT+17)*Y(5)+A(IPT+18)*Y(6) 

IF(NTEMP.GT. 7)  GO  TO  150 
DO  140  1=1,6 
B(I,NTEMP)=Y(I) 

140  B(I,NTEMP+7)=YBAR(I) 

NTEMP=NTEMP+1 
GO  TO  170 


C RESET  STORED  VALUES  OF  Y AND  YBAR  IN  PREPARATION  FOR  NEXT  INTEGRATION 
C STEP. 

150  DO  160  1=1,6 
B(I,1)=B(I,2) 

B(I,2)=B(I,3) 

B(I,3)=B(I,4) 

B(1,4)=B(1,5) 

B(I,5)=B(I,6) 

B(I,6)=B(I,7) 

B(I,7)=Y(I) 

B(I,8)=B(I.9) 

B(I,9)=B(I,10) 

B(I,10)=B(I,11) 

B(I,11)=B(I,12) 

B(I,12)=B(I,13) 

B(I,13)=B(I,14) 

160  B(I,14)=YBAR(I) 

170  CONTINUE 


1 


DIMENSION  B(4,1A) ,D(4,14) ,X(4) ,Y(4) ,XBAR(4) ,YBAR(4) ,XM(4) ,YM(4) , 

1 XMBAR(4) ,YMBAR(4) ,P(4) ,Q(4) ,C(4) ,F(4) ,PMNUSC(4) , 

2 QMNUSF(4) ,A(1430) 

EQUIVALENCE  (X(l) ,X1) 

EQUIVALENCE  (Y(1),Y1) 

EQUIVALENCE  (XM(1),XM1) 

EQUIVALENCE  (YM(1) ,YM1) 

EQUIVALENCE  (XMBARU)  .XMBARl) 

EQUIVALENCE  (YMBAR(l) YMBARl) 

IPT=-12 


C LOOP  OVER  DEPTH  IN  CURRENT-STEP-SIZE  REGION. 

DO  170  N=NSTART,NSTOP 
IF(IT.EQ.4)  GO  TO  115 
IT=NTEMP-4 
115  DO  120  1=1,4 

C SET  PREDICTORS  P(I)  , Q(I)  . 

P(I)=B(I,IT)+COEFFl*(2,0I>4-00*(B(I,IT+10)+B(I,IT+8))-B(I,r'"9)) 

Q(I)=D(I,lT)+COEFFl*(2.0D4-00*(D(I,IT+10)+D(I,lT+8))-D(I,Ii.-9)) 

C SET  MODIFIED  PREDICTORS  XM(I)  , YM(I)  . 

XM(I)=P(I)-.9256198347107438D+00*PMNUSC(I) 

120  YM(I)=Q(I)-.9256198347107438EH-00*QMNUSF(I) 

IPT=IPT+13 

XMBAR(1)=A(1PT)*XM(1)+A(IPT+1)*XM(2)+A(IPT+2)*XM(3) 

Y>fBAR(l)=A(IPT)  *YM(1)+A(IPT+1)*YM(2)+A(IPT+2)*YM(3) 
XMBAR(2)=A(IPT+3)*XM(1)+A(IPT+4)*XM(2)+A(IPT+5)*XM(3) 

1 +A(1PT+6)*XM(4) 

YMBAR ( 2) =A ( IPT+3) *YM ( 1 ) +A ( I PT+4 ) *YM ( 2) +A ( lPT+5 ) *YM ( 3) 

1 +A(IPT+6)*YM(4) 

XMBAR(3)=A(IPT+7)*(XM(3)-XM(1))+A(IPT+8)*XM(4) 
YMBAR(3)=A(1PT+7)*(YM(3)-YM(1))+A(IPT+8)*YM(4) 
XMBAR(4)=A(IPT+9)*XM(1)+A(IPT+10)*XM(2)+A(IPT+11)*XM(3) 

1 +A(IPT+12)*XM(4) 

YMBAR(4)=A(IPT+9)*YM(1)+A(IPT+10)*YM(2)+A(IPT+11)*YM(3) 

1 +A(IPT+12)*YM(4) 

DO  130  1=1,4 

C SET  CORRECTORS  C(I)  , F(I)  . 

C(I)  = .125D+00*(9.0D+00*B(I,IT+3)-B(I,IT+l)-l-COEFF2*(XMBAR(I)+2.0D+0 
1 *B(I,IT+]0)-B(I,IT+9))) 

F(I)=.125D+00*(9.0D+00*D(l,IT+3)-D(I,IT+l)+COEFF2*(YMBAR(I)+2.0D+0 
1 *D(I,IT+10)-D(I,IT+9))) 

PMNUSC(I)=P(I)^C(I) 

QMNUSFU)=Q(I)-F(I) 

C SET  SOLUTION  VECTORS  AT  NTH  DEPTH. 

X(I)=C(I)+.07438016528925620D+00*PMNUSC(I) 

130  YU)=F(I)  + .07438016528925620D+00*QMNUSFU) 

IF(N.EQ.NEND)  RETURN 


C SET  DERIVATIVES  OF  SOLUTION  VECTORS  AT  NTH  DEPTH. 

XBAR(1)=A(IPT)*X(1)+A(1PT+1)*X(2)+A(IPT+2)*X(3) 

YEAR ( 1) =A ( IPT) *Y ( 1) +A ( IPT+1 ) *Y ( 2) +A ( IPT+2) *Y ( 3) 
XBAR(2)=A(1PT+3)*X(1)+A(IPT+A)*X(2)+A(IPT+5)*X(3) 

1 +A(IPT+6)*X(4) 

YBAR(2) =A(IPT+3) *Y (l)+A(IPT+4) *Y (2)+A(IPT+5) *Y (3) 

1 +A(IPT+6)*Y(4) 

XBAR(3)=A(IPT+7)*(X(3)-X(1))+A(IPT+8)*X(A) 

YBARO)=A(IPT+7)*(YO)-Y(lb+A(IPT+8)*Y(4) 

XBAR(4)=A(IPT+9)*X(1)+A(IPT+10)*X(2)+A(IPT+11)*X(3) 

1 +A(IPT+12)*X(4) 

YBAR(4)=A(IPT+9)*Y(1)+A(IPT+10)*Y(2)+A(IPT+11)*Y(3) 

1 +A(IPT+12)*Y(4) 

IF(NTEMP.GT.7)  GO  TO  150 
DO  140  1=1,4 
B(I,NTEMP)=X(I) 

D(I,NTEMP)=Y(I) 

B ( I , NTEMP+7 ) =XBAR ( I ) 

140  D(I,NTEMP+7)=YBAR(I) 

NTEMP=NTEMP+1 
GO  TO  170 

C RESET  STORED  VALUES  OF  X , Y AND  XBAR  , YBAR  IN  PREPARATION  FOR  NEXT 
C INTEGRATION  STEP. 

150  DO  160  1=1,4 
B(I,1)=B(I.2) 

B(I,2)=B(I,3) 

B(I,3)=B(I,4) 

B(I,4)=B(I,5) 

B(I,5)=B(I,6) 

B(I,6)=B(I,7) 

B(I,7)=X(I) 

B(I,8)=B(I,9) 

B(I,9)=B(I,10) 

B(I,10)=B(I,11) 

B(1,11)=B(I,12) 

B(I,12)=B(I,13) 

B(I,13)=B(I,14) 

B(I,14)=XBAR(I) 

D(I,1)=D(I,2) 

D(I,2)=D(I,3) 

D(I,3)=D(I.4) 

D(I,4)=D(I,5) 

D(I,5)=D(I,6) 

D(I,6)=D(I,7) 

D(I,7)=Y(I) 

D(I,8)=D(I,9) 

D(I,9)=D(I,10) 

D(I,10)=D(I,11) 

D(I,11)=D(I,12) 

D(I,12)=D(I,13) 

D(I,13)=D{I,14) 

160  D(1,14)=YBAR(I) 

170  CONTINUE 


K»H1 


C NORMALIZATION. 

AMXINV=1.0D+00/DMAXl(DABS(B(7,l) ) ,DABS(B(7,2) ) , 

1 DABS{B(7,3) ), DABS (B{ 7,4) ) , 

2 DABS(B(7,5) ) ,DABS(B(7,6) ) ) 
DO  165  1=1,6 

PMNUSC (I) =PMNUSC (I) *AMXINV 
Y(I)=Y{I) *AMXINV 
YBAR(I) =YBAR(I) *AMXINV 
DO  165  J=l,14 
165  B(J,I)=B{J,I) *AMXINV 


WAVELENGTH  (km)  PHASE  VELOCITY  PERIOD  (sec) 


MODE  NUMBER 


150  DO  160  1=1,6 

B(1,I)=SNGL(B(2,I)) 
B(2,I)=SNGL(B(3,I)) 
B(3,I)=SNGL(B(4,I)) 
B(4,I)=SNGL(B(5,1)) 
B(5,I)=SNGL(B(6,I)) 
B(6,I)=SNGL(B(7,I)) 

r 

I ‘ B(7,I)=SNGL(Y(I)) 

I B(8,1)=SNGL(B(9,I)) 

B(9,I)=SNGL(B(10,I)) 
B(10,I)=SNGL(B(11,I)) 
B(11,I)=SNGL(B(12,I)) 
B(12,I)=SNGL(B(13,I)) 
B(13,I)=SNGL(B(14,I)) 
160  B(14,I)=SNGL(YBAR(I)) 


max 


PERIOD  ^65  SEC 


GENERATION  OF  COMPLETE  THEORETICAL  SEISMOGRAMS  FOR 


SH  III 

Enzo  Mantovani* 


Publication  No.  1718 

Institute  of  Geophysics  and  Planetary  Physics 
University  of  California,  Los  Angeles 

Publication  No.  000,  authorized  by  the 
Director  of  the  Geodynamics  Project  (CNR) , Italy 

Sponsored  by 

Advanced  Research  Projects  Agency  (DOD) 

ARPA  Order  No.  3291 

Monitored  by  AFOSR  Under  Contract  #F49620-67-C-0038 


4 I 

*Osservatorio  Geofisico,  University  di  Siena,  Italy  I 


RESEARCH  NOTE 


GENERATION  OF  COMPLETE  THEORETICAL  SEISMOGRAMS 
FOR  SH  III 
Enzo  Mantovani 

Siimmary 

Earlier  efforts  to  generate  the  entire  theoretical 
seismograms,  including  both  body  and  surface  waves  for 
realistic  sources  buried  in  a radially  heterogeneous 
anelastic,  spherical  Earth,  are  extended  to  include  the 
summation  of  sixteen  modes . The  comparison  between  a 
real  seismogram  and  theoretical  time  series,  relative  to 
different  attenuation  models  in  the  upper  mantle,  yields 
information  concerning  the  anelasticity  under  the  Pacific 


Ocean. 


Introduction 


Complete  seismograms  which  include  both  body  and  surf^ 
ce  waves  for  a spherical,  anelastic,  radially  heterogeneous 
Earth  have  been  generated  by  simple  inverse  Fourier  treuv- 
sformation  of  the  propagating  fundamental  and  higher-mode 
surface  waves  (Nakanishi  et  al.,197T,  Mantovani  et  al.  ,1977)« 
Efficient  computational  algorithms  permit  us  to  deeG.  with 
highly  realistic  models  of  the  Earth,  whose  radial  hetero- 
geneity is  approximated  by  200  layers,  and  to  employ  2000 
frequency  points  per  radial  mode. 

Theoretical  seismograms 

The  tests  we  report  here  were  carried  out  on  a CIT-tl; 
oceanic  struct\iro  (Anderson  and  Toksoz,  1963).  ^he  first 
sixteen  Love  wave,  or  torsional  modes  were  used,  and  for  e- 
ach,  the  dispersion,  attenuation  and  excitation  are  computed 
down  to  a minimum  period  of  1 second.  The  source  model  we 
have  adopted  is  a dip-slip  displacement  dislocation  on  a ver 
ticeG.  fault  at  a depth  of  T80  km. 

In  Fig,  1 we  illustrate  the  improvement  in  the  details  of 
the  theoretical  seismograms  as  the  number  of  radial  modes 
used  in  the  synthesis  increases.  From  top  to  bottom,  as  the 
n\imber  of  modes  increases, the  amplitude,  clarity  and  short 
period  content  are  seen  to  Increase  for  each  arrival  in  a 
manner  corresponding  to  the  depth  of  penetration  of  the  mo- 
des and  the  depth  at  which  the  arrival  is  produced:  the  first 
six  modes  appear  sufficient  to  reproduce,  almost  completely, 
the  Sa  svirface  wavea  train,  that  is,  the  late -arriving,  domi 
nant  energy  whose  penetration  is  limited  to  approximately 
400km  of  depth  (Kausel  et  al,,1977);  in  the  second  trace, 
synthesized  from  the  first  eleven  modes,  also  the  shallow 
body  wave  phase  SS,  having' a group  velocity  of  about  5.7 
km/sec>  appears;  to  be  almost  completely  formed;  in  the  last 
trace  , synthesized  from  16  modes,  the  deeper  body  wave  ph^ 
ses  S and  SS,  i.e.,  the  first  bursts  on  the  time  series  be- 
gin to  appear  as  two  resolvable  pulses,  but  are  still  devel^ 
ping  in  amplitude  and  short-period  content. 


The  core  reflections  ScS  and  sScS,  whose  penetration  is  much 
greater  than  that  of  the  first  16  modes,  at  the  predominant 
periods  comprising  these  pulses,  are  still  unrecognizable 
in  the  seisqiograms. 

In  Fig*  2 we  show  a range  of  results  out  to  an  epicen_ 
tral  distance  of  8000  km;  the  direct  waives  and  those  reflec- 
ted by  the  surface  of  the  Earth  are  recognizable  in  almost 
every  case;  the  waves  reflected  by  the  core  however  are  re- 
presented only  by  their  long  period  components  ( see  traces 
for  distances  of  3000  and  4000  km,).  Relative  to  the  corre- 
sponding figure  in  Part  II  (Mantovani  et  al.,1977),  it  is 
possible  to  note  a general  iipprovement  in  the  impulsive 
shape  and  resolution  of  individual  body  waves  and  is  parti- 
cularly clear  at  short  distance  in  those  traces  for  the 
15-100  WWSSN  instrument. 

The  comparison  between  experimental  and  theoretical  sei 
smograms  for  30-100  WSSN  instruments  is  given  in  Fig, 3* 

Details  concerning  the  event  and  a preliminary  compari- 
son between  theoretical  and  experimental  time  series  using 
surface  waves  is  reported  by  Kausel  et  al,  (1977),  Our  pur- 
pose here  is  to  continue  the  discussion,  begun  in  Part  II, on 
the  qualitative  aispects  of  the  theoreticeO.  - experlmentaJ. 
comparison  as  more  and  more  higher  modes  are  included  in  the 
synthesis  of  the  theoretical  traces.  We  also  report  on  some 
tests  of  the  sensitivity  of  the  amplitudes  and  arrival  times 
of  the  pulse  bursts  to  the  model  of  the  Earth’s  intrin 

sic  attenuation  and  to  the  source  parameters. 

In  Part  II  the  variation  of  relative  amplitudes  of  body 
and  surface  waves  caused  by  chauiges  in  crustal  emd  upper-man 
tie  attenuation  was  noted,  and  it  wais  suggested  that  this 
effect  might  be  a useful  tool  in  the  investigation  of  the 
Earth's  Intrinsic  anelasticity.  We  have  used  two  models  of 
attenuation.  Model  I is  based  on  model  MM8  of  Anderson  et 
al,(l965).  This  model  wan  obtained  from  measurements  of  the 
decay  of  Rayleigh  and  Love  waves  traveling  globe  circling 
paths;  thus  it  cannot  be  expected  to  be  an  accurate  represen- 
tation of  the  anelasticity  in  a specific  oceanic  region. 


Model  II  is  based  on  that  given  by  Mitchell  (1976)  derived 
from  inversion  of  Raylei^-wave  attenuation  measurements  in 
the  "older  portions  of  the  Pacific  crust  and  upper  mantle". 
Models  I and  II  are  shown  in  Pig. 4. 

Since  the  average  anelasticity  in  the  first  400  km  is  very 
similar  for  modni  I ( Q = 0.00825  ) and 


model  II  ( Q = 0.00774 


) the  Sa  waves 


have  about  the  same  attenuation.  The  ratio  of  amplitudes  of 
Sa/S  is  about  the  same  for  both  models  since  amplitudes  of 
S are  only  slightly  aiffected  by  variations  of  anelasticity 
in  a small  portion  of  their  path.  Since  Sa  waves  sample  the 
uppermost  400  km  of  the  structvire  almost  uniformly,  they  are 
more  sensitive  to  the  average  attenuation  in  this  region  thain 
to  its  detailed  depth  distribution. 

A comparison  of  traces  A and  B (Fig. 3),  with  the  help  of  TA- 
BLE 1,' indicates  that,  for  a source  at  180  km,  the. relative 
amplitudes  of  S and  Sa  are  not  so  different  for  model  I and 
II.  The  striking  featvire  of  trace  B,  relative  to  A,  is  the 
remarkable  increase  in  the  amplitudes  of  SS  arrivals.  It  is 
reasonable  to  suppose  that  this  effect  is  connected  with  the 
lower  attenuation  of  model  II  between  about  250  and  700  km, 
that  is  in  that  part  of  the  structure  in  which  the  SS  rays 
have  a significant  portion  of  their  path. 

A comparison  of  the  well  pronounced  SS  phase  in  trace  B,  with 
even  the  largest  amplitudes  around  the  saijie  arrival  in  the 
experimental  series,  leads  us  to  suspect  that  model  II  has 
a drop  in  attenuation  which  is  too  large  between  2 50  and  700 
km  of  depth  (Fig.4). 


On  the  basis  of  the  results  shown  in  Fig,  1,  we  infer  that  the 
above  observations  are  little  affected  by  the.  finite  nvimber 
of  radial  modes  used  in  this  study  since  the  phase  SS  appears 
to  be  converging  both  to  a definitive  wavelet  form  and  ampli- 
tude. 

A final  note  in  this  series,  which  will  describe  theoretical 
seismograms  generated  by  even  more  radial  modes,  will  give 
more  information  on  this  point. 

To  Investigate  further  the  effect  produced  by  changes  in  the 
average  anelastlcity  of  the  crust  and  upper  mantle,  we  gene- 
rated time  series  (Fig. 3, trace  C)  corresponding  to  an  anela- 
Bticity  distribution  derived  from  model  I,  but  with  the  pha- 
se attenuation  Bg  , reduced  to  0,000700  sec/km  from  the  top 
of  the  low-velocity  channel  down  to  the  400  km  discontinuity. 
Comparison  of  traces  C and  A shows  an  appreciable  decrease 
in  the  amplitude  ratio  S/Sa,  and  a significant  increase  of 
the  energy  in  the  tail  of  the  Sa  wave  train. 

As  a limiting  case  for  the  anelastlcity  tests  we  have  also 
exhibited  the  time  series  (Pig.  3,  trace  D)  for  a perfectly 
elastic  Earth,  In  this  case  ,for  reasons  reported  in  Part  II, 
the  surface  waves  completely  dominate  those  body  waves  that 
penetrate  deeply. 

The  phases  SS  show  a less  drastic  decrease  in  amplitude  be- 
cause they  are  relatively  strongly  affected  by  the  high  atte- 
nuation in  the  upper  mantle. 

To  give  a representation  of  the  effect  , on  the  time  series, 
produced  by  variations  in  source  depth,  we  present  theoreti- 
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res  as  traces  but  with  a loeal  depth  of  140  km* 

We  see  that  variation  of  the  source  depth  has  an  effect 
on  the  arrival  times  and  the  ao^litudes  of  body  waves;  the 
feature  of  the  time  series .which  is  most  sensitive  to  the 
eo\irce  depth,  at  leaist  for  the  focal  and  structural  parame 
ters  we  used,  is  the  amplitude  of  the  SS  phase. 
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FIGURE  CAPTIONS 


Representation  of  the  effect  of  adding  more  higher  (radial) 
modes  in  the  synthesis  of  theoretical  sismograms  recorded 
at  the  surface  of  the  Earth,  Source  mechanism  is  a dip 
slip  displacement  dislocation  on  a vertical  fault  at  a depth 
of  180  km.  Epicenter-station  separation  is  6788  km.  The 
model  we  have  used  for  the  intrinsic  anelasticity  of  the 
Earth  is  an  adaptation  (Fig, 4)  of  model  MM8  (Anderson,  Ben 
Menahem  & Archambeau  , 1965)  to  our  oceanic  structure. 

Torsional-wave  response,  at  the  surface  of  the  Earth,  to 
a dip-slip  displacement  dislocation  on  a vertical  fault, 
at  a depth  of  180  km.  The  expected  body-wave  arrival  times 
are  indicated  for  S(  •) , SS  (A),  scS  (▼),  sS  (■),  sScS 
(♦)  , See  caption  of  Pig.  1 for  anelasticity  model. 

Comparison  of  theoretical  and  experimental  seismograms 
for  30-100  WWSSN  instruments,  6788  km  from  the  epicenter; 
short  traces  are  copies  of  the  long-period  experimental 
record  (SH)  measured  at  Honiara  (hnr)  for  the  event  occur- 
ring at  the  foot  of  the  Kamchatka  peninsula  at  14:30:30.3 
GCT  on  1964  December  26  (see  Kausel  et  al,,  1977  for  full 
details  on  this  event,  our  theoretical  treatment  of  the 
source  mechanism,  and  the  analysis  of  the  Sa  portion  of 
the  time  series).  Traces  denoted  by  capital  letters  are 
theoretical  seismograms;  trace  A represents  the  anelasti 
city  Model  I ( Fig,  4 ) ; trace  B,  anelasticity  Model  II 
(Fig,  4);  trace  C,  MD.del  I with  the  phase  attenuation  Bg 
reduced  to  0,000700  sec/km  from  the  top  of  the  low  velo- 
city channel  down  to  the  ''400-km''  discontinuity;  trace  D 
does  not  include  attenn^ticn  Sw 
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E,  F,  6,  and  H represent  the  same  models  respectively  as 
A,  B,  C and  D but  with  the  source  depth  reduced  from  180 
to  140  km. 

Fig,  4 S-wave  velocity  structure,  transverse-wave  phase  attenua- 
tion Bg,  and  intrinsic  attenuation  0*”\  obtained  from  the 
adaptation,  to  our  velocity  structure,  of  two  anelasticity 
models  reported  in  the  literature;  (I)  world-wide  average 
model  (Anderson,  Ben  Menahcm,  Archambeau,  1965),  (ll)  0- 
ceanic  model  (Mitchell ,1976) , Model  I is  given  by  solid 
lines.  Model  II  by  dotted  lines  . 


TABLE  1 

Amplitude-ratios  S/S^  and  SS/S  measured,  peak  to  peak,  on  the 

a di 

traces  A,  B,  C,  E,  P and  6 in  Pig,  3*. 

180  km  140  km  ■ 

ABC  E P G 


s/s^ 

0.66 

0.62 

0.54 

0.59 

0.50 

0.43 

ss/s^ 

0.27 

0.39 

0.25 

0.48 

0.58 

0.41 
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DEFINITIONS  OF  QUANTITIES  USED  IN  THIS  PAPER 


G 

u;-. 
uQz 
' t, 

ei 

Va, 

(T, 

01 

cr-^ 

H 

K 

A 

K 
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Angular  deflection  of  seismometer  boom  in  radians 
Angular  deflection  of  galvanometer  mirror  in  radians 

Angular  frequency  of  seismometer  free  period,  =Xr^f~^  , (radians/sec) 
Angular  frequency  of  galvanometer  free  period,  , (radians/sec) 

Seismometer  damping  parameter,  = U.\ 

Galvanometer  damping  parameter,  = VN2.«tC2^ 

Seismometer  damping  constant  (1.0  = critical  damping) 

Galvanometer  damping  constant  (1.0  = critical  damping) 

Seismometer  coupling  constant 
Galvanometer  coupling  constant 

= 'T'/'jI  , coupling  constant  between  seismometer  and  galvanometer 
(1.0  = direct  coupled) 

Seismometer  boom  mass  (kilograms) 

Distance  between  center  of  gravity  of  mass  and  rotation  point 
of  seismometer  boom  (meters) 

Moment  of  inertia  of  seismometer  boom  (kg-m"*) 

Displacement  of  the  ground  (meters) 

= G;JK.  , magnitude  of  the  force  impressed  on  the  boom  by 
*the‘ cal ibration  current  (Newtons) 

Calibration  current  (mi  1 1 i amperes ) 

Motor  constant  of  the  calibration  coil  (Newtons/Ampere) 

Height  of  the  impulse  response  as  measured  on  the  record  (milli- 
meters) 


VHi ^ Heaviside  step  function  of  time 
S(,t)  Dirac  delta  function  of  time 

TAi/  Pseudoparameter  forms  of  the  seismometer  and  galvanome- 

' ter  free  periods  and  dampings 


The  time  at  which  the  calibration  current  was  switched  on 

An  apparent  amplitude  factor  that  is  determined  in  the  inversion 

A constant  term  to  be  subtracted  from  the  digitized  data  to 
best  agreement  with  theoretical  impulse  responses 

A linear  trend  to  be  subtracted  from  the  digitized  data  to 
give  best  agreement  with  theoretical  impulse  responses 

A rotation  to  be  applied'  to  the  impulse  response  to  correct 
for  incorrect  digitizing  of  the  impulse 
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IMPULSE  RESPONSE  INVERSION 


1 


Quantitative  time  series  analysis  using  the  World-Wide  seismic 
network  records  (WWSSN)  requires  that  corrections  be  made  for  the 
response  of  the  seismograph.  In  the  usual  method  of  calculation  of 
instrumental  corrections  the  seismograph  Is  considered  as  a black  box. 
From  Fourier  analysis  of  the  response  to  a known  transient  impulse 
(usually  a step  function  in  acceleration  applied  to  the  seismometer 
mass)  corrections  can  be  obtained.  In  general,  at  frequencies  in  the 
neighborhood  of  the  peak  instrumental  response  this  is  an  accurate 
method  of  correction.  However,  periods  shorter  than  the  seismometer 
free  period  are  not  strongly  excited  in  the  impulse  response  and  sig- 
nificant errors  in  the  computation  of  ground  motion  may  occur. 

Studies  involving  short  period  crustal  and  higher  mode  surface 
waves  have  obliged  us  to  seek  an  alternative  method  of  determining  these 
corrections,  In  order  to  reduce  errors  caused  by  Instrumental  uncer- 
tainties. The  parameters  of  the  seismometer  system  have  been  determined 
directly  using  a linear  inversion  method.  Mitchell  and  Landlsman's 
(1969)  analytical  expressions  can  then  be  used  to  derive  the  phase  and 
amplitude  corrections  at  all  periods.  This  method  uses  a stabilized 
1 eas t- squa res  technique  which  allows  more  parameters  to  be  used  while 
maintaining  convergence  in  the  inversion  process.  Arbitrary  values  of 
the  coupling  constant  and  the  seismometer  period  can  be  used  in  this 
form  of  the  calculation.  Also,  the  scale  factor  determined  in  the  in- 
version can  be  used  to  obtain  an  expression  for  seismograph  magnifica- 
tion that  is  independent  of  the  rated  value. 
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The  procedure  for  inversion  ( Jackson , 1 9 72 ) involves  solving  a sys- 
tem of  linear  equations  for  the  corrections  'to  the  model  para- 

meters. In  our  case  these  parameters  are  the  seismograph  instrumental 
constants.  The  equations  are  set  equal  to  the  data  residuals 
which  are  the  residuals  between  the  data  and  the  predictions  from 
the  model.  The  linear  equations  take  the  form 


itrix  partial  derivatives  of  the 

data  residuals  with  respect  to  the  model  parameters.  The  I east- squares 


where  the  mal 


solution  results  from  the  straightforward  minimization  of  the  error 
residuals  (Mitchell  and  Landisman,  1969)-  They  noted  that  the  in- 
clusion of  too  many  parameters  tended  to  cause  instabilities  in  the 
inversion.  Thus  certain  parameters,  such  as  the  scale  factor  and 
impulse  'on'  time,  had  to  be  empirically  adjusted  in  the  inversion. 

In  our  experience  with  digitizing  actual  impulse  responses  there  is 
error  introduced  due  to  baseline  choice,  and  so  a constant,  linear 
trend,  and  rotation  parameters  were  introduced.  Although  all  the 
partial  derivatives  have  linearly  Independent  components,  in  the 
presence  of  noise  not  all  of  them  can  be  resolved,  and  an  unstable 
^inversion  resulted  using  1 eas t- squares  method  with  the  nine  parameters 
listed  in  figure  1.  The  inversion  was  stabilized  by  decreasing  the 
degrees  of  freedom  in  the  inversion  (Wiggins,  1972).  The  matrix 

of  partial  derivatives  is  first  normalized  with  respect  to  the  data 
standard  deviations  and  a priori  uncertainties  assigned  to  the  model 
parameters.  The  matrix,  now  in  normal  form  (Jackson,  1972),  is  factored 
according  to  the  Lanczos  factorization  (Lanczos  1961)  Into  eigenvalue 
and  eigenvector  matrices.  In  this  form  eigenvalues  less  than  1.0 
correspond  in  some  sense  to  parameter  combinations  that  cannot  be 
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resolved  from  the  data,  and  so  these  eigenvalues  are  set  to  zero, 
removing  these  degrees  of  greedom  from  the  inversion.  In  general 
the  number  of  independent  degrees  of  freedom  in  the  Inversion 
when  this  was  done  was  six  or  seven.  All  the  parameters  necessary 
to  obtain  phase  and  magnification  Information  were  obtainable, 
but  careful  digitization  was  necessary  to  insure  good  results. 

COMPUTATIONS 

The  actual  computation  of  the  partial  derivatives  was  done 
using  an  analytical  formula  for  the  inpulse  response  as  a func- 
tion of  the  parameters.  Mitchell  and  Land i sman( 1 969)  used  the 
inverse  Fourier  transform  of  the  frequency  domain  representation 
of  the  impulse  response  to  calculate  the  partial  derivatives. 

While  this  has  certain  advantages,  I felt  that  the  number  of  Fou- 
rier transforms  needed  per  inversion  would  be  a significant  cost 
factor  even  with  FFT,  and  so  a method  was  followed  as  illustrated 
in  figure  2.  Starting  with  the  representation  of  the  seismometer 
and  galvanometer  as  coupled  damped  harmonic  oscillators,  the 
equation  of  motion  of  the  rotation  of  the  galvanometer  Is  solved 
by  inverse  Laplace  transform.  The  magnitude  of  the  impulse  given  to 
the  mass,  K , is  just  . Jarosch  and  Curtis  (1973)  noted 

that  the  denominator  cannot  be  factored  unless  the  term  containing 
(pZ.  is  set  to  zero.  When  this  is  done  the  pseudo- parameter 
equation  can  be  solved  uniquely  to  obtain  new  values  of  w',  , 

and  we  call  the  new  values  of  the  'pseudoparameters'. 

As  can  be  seen  from  the  figure,  the  effect  of  this  conversion 
is  small  for  small,  but  may  be  significant  for  larger  . 
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Furthermore,  the  pseudoparameters  become  complex  for  larger 

We  can  obtain  finally  an  analytic  expression  for  the  impulse  response: 


^ 

whJre-  |<"T|^7^rr'  ^ 

This  will  be  in  general  comp  1 ex- va 1 ued  for  complex  parameters; 

' the  real  part  corresponds  to  the  impulse  response  function. 

If  the  seismometer  or  ^Ivanometer  becomes  cri  t i cal  1 y damped, 
we  must  remember  that  Ki-O  or  Ka.  = 0 , and  so 

-i  ; CC5  1 ^ cos  Kai  =i  • 

For  the  case  when  the  seismometer  or  galvanometer  is  overdamped, 

K|  Ki  becomes  complex,  and  so  the  appropriate  sin  and 

cos  are  replaced  by  sinh  and  cosh.  In  figure  3 are  some  impulse 
responses  calculated  using  the  above  formula.  The  reference  impulse 
response  corresponds  in  some  way  to  a standard  WWSSN  instrument. 

Some  notion  of  what  the  partial  aerivatives  look  like  can  be 
gained  from  this  figure.  The  effect  of  changing  the  seismometer 
free  period  can  be  seen  to  be  very  similar  to  the  effect  of  chang- 
ing the  gain,  or  amplitude  factor  of  the  impulse  response.  This 
effect  can  be  seen  in  the  Lookout  Mountain  calibration  tests,  in 
which  two  different  impulse  responses  from  the  instrument  were 
. inverted  independently.  The  seismometer  free  period  values  show 
somev;hat  larger  discrepancies  than  the  other  parameters.  The  theo- 
retical values  for  the  parameters  compare  very  well  with  those 
determined  experimentally  using  the  usual  tests  for  damping  and 


z:\ 


M 


-5- 


free  period.  In  the  case  of  the  galvanometer  constants,  there  is 
a very  good  argument  for  determining  the  constants  through  inversion 
rather  than  by  other  methods,  since  this  seems  to  be  at  least 
as  accurate  and  does  not  disturb  the  instrument  in  any  way. 

The  seismometer  constants,  especially  the  free  period,  are  probab- 
ly better  determined  by  other  experimental  metliods,  and  constrained 
in  the  inversion.  Once  the  instrumental  constants  are  determined, 
the  instrumental  effects  on  amplitude  and  phase  can  be  determined 
theoretically.  Following  Mitchell  and  Landisman  (.1969)  the  phase 


can  be  wr i tten  as : 


vtO;*  t-  10^-  ♦-  to 


I W I 


The  convention  here  is  that  'S  positive,  which  corres- 

ponds to  positive  instrumental  group  delay  time.  The  magnification 
of  the  instrument  can  similarly  be  expressed  in  terms  of  the  para- 


meters of  the  instrument: 


I'p) 


is  the  inversion  gain  constant,  tp  is  the  calibration  coil 
current,  and  Q is  the  motor  constant  of  the  coil.  For  WWSSN 
stations  there  is  an  empirical  formula  for  the  gain  at  the  seis- 
mometer free  period: 


Vir 


for  T,  = 1 5 seconds  ; 


Gif 


for  =30  seconds 


IS  the  height  of  the  impulse  response  in  millimeters.  This 


formula  was  checked  against  the  analytical  formula  for  gain  and 
agreement  to  a few  percent  or  better  was  found.  However,  the 
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empirical  formula  assumes  all  the  constants  are  at  the  rated  values  ,, 

and  so  a badly  tuned  instrument  would  give  erroneous  gains. 

A long  period  WWSSN  station,  NDI-LPZ,  was  used  in  inversion  to 
get  instrumental  parameters.  The  amplitudes  in  the  impulse  response 
must  be  multiplied  by  UJ"^  ^ to  get  the  ground  magnification.  This 
means  that  the  gains  at  the  long  period  end  are  well  determined, 
where  the  original  impulse  had  a strong  signal.  However,  below 
30  seconds  period  the  amplitudes  become  relatively  worthless  and 
it  is  advisable  to  use  the  theoretical  values.  The  same  consider- 
ations hold  for  the  phases.  An  additional  source  of  error  in  the 
phase  of  the  impulse  is  introduced  because  the  'on'  time  of  the 
impulse  is  not  known  a priori , but  must  be  estimated  from  the  re- 
cord. An  error  of  one  second  or  more  is  common  here.  This  leads 
to  a systematic  bias  in  the  phases  determined  as  can  be  seen  in 
figure  5,  where  the  empirically  determined  'on'  time  is  wrong 
by  one  second.  No  bias  is  evident  in  the  phases  derived  using 
the  'on'  time  from  inversion. 


A note  concerning  baseline  determination  is  relevant  here. 
Although  it  was  not  possible  in  the  inversion  to  use  the  rotation 
as  an  independent  parameter,  the  rotation  is  quite  important,  es- 
pecially for  very  large,  steep  impulses,  and  a wrong  baseline  will 
distort  the  constants  obtained  and  resultant  phases  quite  badly. 
However,  dravying  the  correct  baseline  is  not  just  a matter  of  con- 
necting the  ends  of  the  trace  and  making  sure  the  impulse  does  not 
occur  in  a disturbed  part  of  the  record.  Figure  7 shows  an  impulse 
on  our  instrument.  The  minute  marks  line  up,  which  makes  it  easy 
to  notice  the  offset  of  the  minute  mark  near  the  top  of  the  pulse. 

gnment.  This  effect  may  be 
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common  among  WWSSN  Instruments,  as  noted  by  James  and  Linde  (1971). 
What  I did  for  this  was  to  rotate  the  impulse  on  the  digitizer  un- 
til the  minute  marks  on  the  pulse  fell  at  equal  intervals  on  the 
digitizer  scale.  This  meant  that  the  digitizer  scale  was  then 
aligned  parallel  to  the  direction  of  motion  of  the  galvanometer 
light  beam.  Digitizing  the  impulse  in  this  manner  gave  satisfactory 
results  even  for  large  impulse  responses. 

Some  conclusions  I have  reached  in  this  work  are  that  the  appli- 
cation of  a generalized  inversion  scheme  to  instrument  impulse 
responses  allows  in  inclusion  of  more  parameters  without  sacrifi- 
cing the  stability  and  validity  of  the  inversion.  In  particular, 
the  use  of  the  amplitude  factor  allows  an  independent  measure  of  the 
magnification  to  be  macie.  The  Inclusion  of  the  'on'  time  means 
that  this  source  of  phase  error  can  be  practiv.ally  eliminated. 

The  constant  and  linear  trend  terms  are  important  because  of  digi- 
tizing considerations  in  which  a baseline  must  be  assumed  and  later 
corrected  for.  Careful  digitizing  is  necessary  to  assure  that  the 
correct  baseline  is  used. 


Noise  on  the  record  during  the  Impulse  response  did  not  seem 
to  be  a significant  factor  in  the  inversion.  In  general  I chose 
only  impulses  that  occured  during  clear,  quiet  times  and  without 
visible  noise  except  for  microseisms.  Examining  the  residuals  that 
remain  after  inversion,  I found  that  much  of  the  error  seemed  to  be 
associated  with  digitizing  errors  on  the  steep  parts  of  the  Impulse. 
This  is  unavoidable  because  of  the  steepness  of  the  pulse  and  the 
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for  good  impulses  that  had  successful  inversions.  These 
errors  did  not  appear  to  cause  much  error  in  the  parameters 
determined  as  shown  by  the  good  agreement  between  indepen- 
dent impulse  inversions  from  the  same  instrument. 

ERROR  ANALYSIS  AND  CONCLUSIONS 

Errors  in  the  determination  of  parameter  values  are  of 
several  different  types.  One  type  of  error  arises  because 
of  the  finite  line  thickness  and  other  digitizing  errors,  and 
because  of  microseisms  and  other  types  of  noise  present  on 
the  record.  These  effects  can  all  be  included  as  noise  in 
the  impulse  response  signal.  Figure  8 shows  the  results  of 
inversion  of  theoretical  and  actual  impulse  responses. 

Random  noise  (white  noise)  was  added  to  one  theoretical 
impulse  and  inversion  was  done.  Convergence  of  error  to  the 
noise  level  was  quite  rapid,  and  after  that  essentially  no 
change  took  place.  Similar  results  were  obtained  for  two 
actual  impulses,  LMO-2  and  NDI-LPZ.  The  digitizing  noise 
level  was  obtained  by  an  estimation  of  the  random  error  to 
be  about  .14  mm  on  the  record  on  the  average.  However,  in 
the  steep  parts  of  the  impulse  the  errors  became  greater,  and 
residuals  of  amplitude  equal  to  about  1%  of  the  impulse 
response  amplitude  (about  1 mm)  were  occasionally  seen. 

Since  the  steep  regions  of  the  impulse  response  are  most 
critical  in  determination  of  the  instrumental  parameters, 
this  upper  bound  was  taken  as  the  representative  noise  level 
In  the  analysis  of  propagation  of  errors.  The  effect  of 
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I 

random  errors  in  the  data  on  the  parameters  was  estimated 

by  computing  the  covariances  of  the  model  parameters  by  the  i 

method  of  Jackson  (I972)C  ] 

The  matrix  \\  is  ~ \J  l\  Vi  , the  general  inverse 
of  the  matrix  of  partial  derivatives.  The  results  obtained 
using  1 mm  as  the  random  error  are  shown  in  table  1.  The  in- 
strumental parameters  for  LMO- 1 and  LMO-2  appear  to  agree 
within  the  uncertainties  given  for  the  parameters,  although 
the  errors  in  the  data  may  not  actually  be  random  as  assumed. 

The  errors  on  the  seismometer  and  galvanometer  free  periods  are 
about  equal  for  the  severely  underdamped  case,  LMO- 1 and  LMO-2, 
but  change  for  the  overdamped  case,  LMO-NEW.  The  impulse  response 
was  activated  by  an  automatic  relay  at  exactly  12-hour  intervals. 

There  appears  to  be  some  difference  between  the  'on'  times  of  the 
impulses  LMO- 1 and  LMO-2  on  the  seismogram,  as  judged  by  measure- 
ments relative  to  the  minute  marks.  This  may  account  for  the 
difference  between  the  Tc>  times  of  LMO- 1 and  LMO-2  impulse 
responses.  0 , the  rotation,  was  constrained  to  equal 

zero  throughout  the  inversions. 

Another  type  of  error  is  not  random,  and  results  because  the 
seismometer  free  period  and  the  amplitude  scale  parameter  p\  can- 


n 1 

• , -10- 

« 

! constraining  in  the  inversion,  since  this  parameter  can 

be  accurately  measured.  Table  2 shows  the  result  of  constraining 
Xj  by  adding  in  the  parameter  eigenvector  to  the  solution. 

Xt  and  have  both  been  made  consistent  by  this  procedure. 

Some  conclusions  that  can  be  drawn  from  this  work  are  that  ap- 
plication of  generalized  linear  inverion  to  seismograph  impulse 
responses  to  derive  the  instrument  constants  yields  useful  and 
consistent  results,  and  makes  possible  more  accurate  estimates 
of  gain  and  phase  shift  at  short  periods.  An  independent  means 
of  deriving  the  gain  as  a function  of  period  is  possible  using 
the  inversion  amplitude  parameter.  Errors  in  phase  due  to  un- 
certainty in  impulse  'on'  time  are  reduced.  There  is  a trade- 
off of  \»  , the  seismometer  free  period  and  f\  , the  gain 

constant,  but  these  can  be  resolved  by  constraining  X > which 
is  in  general  well-known.  Noise  caused  by  digitizing  errors  or 
microseisms  does  not  appear  to  be  a problem  in  this  work.  The 
residuals  were  usually  of  amplitude  less  than  1%  of  the  original 
impulse  amplitude.  It  v;as  necessary  to  take  great  care  in  choos- 
ing the  correct  baseline  for  digitizing.  Application  of  the 
'nversion  technique  to  actual  impulse  responses  yielded  consis-  ' 

ent  results  and  enabled  amplitude  and  phase  shifts  to  be  cal- 
culated at  all  periods. 
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APPENDIX  A 

ANALYTICAL  METHOD  FOR  OBTAINING  PSEUDOPARAMETERS 


From  figure  2 we  have  the  equality  involving  original  and  i 

pseudoparameters  that  must  be  satisfied:  1 

{ S - + ae ( &■•' • ia .s u;. H iff = - ( sM'c lo; 5 H-:  I 

Following  Jarosch  and  Curtis  (1971)  we  equate  powers  of  S 


(1) 

(2) 

(3) 

(4) 

To 

simplify,  let 

(5) 

W - CK. 

(6) 

(7) 

(8) 

t i 

(9) 

We 

can  then  derive 

the  characteristic  equation  in  CL 

(10)  -roy  -<\y)  0?  - |/vr^  = l 


Rather  than  solve  this  equation  numerically  as  Jarosch  and  Cur- 
tis did,  I needed  the  complex  solutions  to  the  equation.  The 
equation  can  be  factored  as 


We  must  solve  ■ Since  there  is  symmetry  between  x,l| 

they  may  be  interchanged: 


(12) 


- lAtS ) (ft-G)  ff'  - PF' 
■'  •2,  ^ Z ’ X Z 


where  y 

(13)  ft:  t 

(14)  G - - IwP-,/(2'7')'''-3 

and 

(15)  ^ - /:3  -''"^3 

(16)  V)  r - j^^~/ 3 - 

(17)  ok'  ^ \r^->  ^ 

(18) 

The  solution  for  ICj*”  is  

( 19)  ' ~ 'j')"- 

* (-E  ±- 

% Z- 

I always  used  the  positive  sign  in  the  root.  If  one  root  was 
the  other  two  corresponded  tou>^andu>\.  It  was  always  possible 
to  distinguish  which  was  which. 


FIGURE  CAPTIONS 


Figure  1:  Schematic  model  of  the  impulse  response  inversion 
program.  The  periods  and  dampings  are  in  pseudoparameter  form 
as  used  in  the  inversion.  Up  to  five  iterations  were  used  in 
convergence  to  the  stabilized  1 eas t- squa res  solution. 

Figure  2;  Equations  for  the  seismometer  and  galvanometer  motions, 
following  Jarosch  and  Curtis  (1973)«  The  effect  of  the  pseudo- 
parameter conversion  is  small  for  small,  but  significant 

for  larger 

Figure  3:  Effect  of  various  parameters  on  the  impulse  response 
function.  The  solid  line  is  the  for  all  figures.  Changing  the 
seismometer  free  period  has  a large  effect  on  the  impulse  am- 
pl i tude. 

Figure  4:  Results  of  inversion  on  our  long  period  vertical  in- 
strument. The  parameters  were  measured  in  s i tu  using  overshoot 
ratio  tests  and  timing  the  free  periods  ¥y  stopwatch. 

Figure  5:  Phase  determination  of  NDI-LPZ  impulse  response,  using 
direct  Fourier  transform  method  and  the  inversion  method.  The 
Fourier  method  is  degraded  at  periods  below  30  seconds.  The 
'on'  time  had  been  picked  by  eye  at  4 seconds  after  the  minute 
mark;  the  inversion  showed  It  actually  occured  at  3 seconds  after 
the  minute.  The  large  black  dots  correspond  to  the  best  esti- 
mate of  the  phase. 

Figure  6:  Amplitude  determination  of  NDI-LPZ  using  the  Fourier 
method  and  inversion  method.  The  curve  on  the  left  corresponds  to 
impulse  amplitudes.  These  become  very  small  below  30  seconds, 
which  accounts  for  the  degradation  of  the  phases  and  amplitudes. 

To  obtain  ground  displacement  magnification,  this  curve  is  mul- 
ti pi  i ed  by 

Figure  7:  Impulse  response  on  the  Lookout  Mountain  (LMO)  LPZ 
record.  There  is  a slight  misalignment  of  the  light  beam  rela- 
tive to  the  drum  as  shov./n  by  the  offset  of  the  minute  mark  on 
the  impulse.  This  minute  mark  was  used  to  help  align  the  impulse 
for  correct  digitization. 

Figure  8:  Convergence  of  the  linear  inversion  method  for  theo- 
retical and  actual  impulse  responses.  For  all  successful  solu- 
tions the  method  converged  within  5 iterations. 
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COMPUTATION  OF  COMPLETE  THEORETICAL 
SEISMOGRAMS  FOR  TORSIONAL  WAVES 

By  a.  H.  Liao,  F.  Schwab,  and  E.  Mantovani 
Introduction 

In  the  study  of  earth  structure,  and  earthquake  source  mechanism,  the  handiest 
and  probably  most  important  information  is  that  provided  by  the  seismic  waves 
recorded  on  seismograms.  Therefore,  the  direct  comparison  between  theorerical 
and  experimental  seismograms  is  a very  appealing  subject  for  most  seismologists. 

The  computation  of  theoretical  seismograms  has  been  studied  by  various  authors. 
Some  tried  to  obtain  long-period  seismograms  by  combining  a few  normal  modes; 
some  tried  to  obtain  a few  body-wave  phases.  The  results  of  the  investigations 
that  led  to  the  present  letter  are  reported  by  Schwab  (1970);  Schwab  and  Knopoff 
(1970,  1971,  1972,  1973);  Kausel  and  Schwab  (1973);  Schwab  and  Kausel  (1976); 
Knopoff  et  al.  (1973);  Knopoff  ef  al.  (1974);  Schwab  et  al.  (1974);  Nakanishi  et  al. 
(1976);  Kausel  et  al.  (1977);  Nakanishi  et  al.  (1976);  Mantovani  et  al.  (1976); 
Mantovani  (1977a);  and  Mantovani  in  preparation).  The  results  of  our  most 
recent  work  now  allow  us  to  generate  complete  theoretical  seismograms  for  torsional- 
wave  motion.  By  saying  complete,  we  mean  that  all  modes  that  exist  for  periods 
above  10  sec  are  included  in  our  seismograms,  i.e.,  that  all  amplitude  and  phase 
information  down  to  a period  of  10  sec  is  included.  This  means  that  the  body-wave 
phases,  as  well  as  surface-wave  arrivals,  are  obtained  with  this  method. 

The  computational  technique  is  somewhat  involved;  however,  anyone  who  is 
interested  in  the  details  can  obtain  program  copies  horn  the  authors.  Briefly,  we 
compute  all  of  the  required  frequency-domain  information  such  as  phase  and  group 
velocity  dispersion,  attenuation,  amplitude  excitation,  apparent  initial  phase,  depth 
of  penetration  (and  if  desired,  partial  derivatives  with  respect  to  structural  param- 
eters) from  the  longest  existing  period  for  each  mode,  down  to  10  sec.  This  requires 
between  90  and  100  radial  modes,  and  is  accomplished  in  a single,  relatively  short 
computer  run.  This  computation  corresponds  to  about  50,000  free  oscillating  modes. 
Theoretical  seismograms  in  the  time  domain  are  then  obtained  by  inverse  Fourier 
transformation  as  described  by  Kausel  and  Schwab  (1973)  and  Schwab  and  Kausel 
(1976).  The  cost  of  computing  the  volume  of  frequency-domain  information,  which 
is  required  to  construct  the  synthetic  seismograms  down  to  a period  as  short  as  10 
sec,  has  been  prohibitively  high  up  to  the  present  time;  but,  with  our  current 
optimization  we  have  succeeded  in  reducing  the  computation  time  to  about  6 min 
on  our  IBM  360/91  computer  (an  expenditure  of  about  $50)  for  a given  earth 
structure.  The  corresponding  time  on  a CDC  7600  computer  is  about  2 min.  [For 
corresponding  times  on  other  computers  see  Schwab  et  al.  (1977)  and  Porter  et  al. 
(in  preparation).]  These  time  estimates  should  be  considered  only  as  upper  bounds; 
relatively  simple  improvements — which  a little  practical  experience  will  be  needed 
to  justify — are  expected  to  cut  these  times  at  least  in  half. 

The  purposes  of  the  present  letter  are  (1)  to  report  that  the  capability  now  e.tbts 
for  computing  realistic,  torsional-wave  seismograms  that  contain  all  of  the  seismic 
energy  and  arrivals  down  to  a period  of  10  sec;  (2)  to  announce  the  availability,  fur 
general  distribution,  of  the  associated  program;  and  (3)  to  exhibit  the  results  of  our 
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initial  testa  of  direct  compariaon  of  these  complete  theoretical  seismograms  with 
an  experimental  record  from  a standard,  long-period  instrument  of  the  WWSSN. 

Examples 

For  our  comparison  of  theoretical  and  experimental  seismograms,  we  have  chosen 
an  event  with  its  epicenter  near  the  foot  of  the  Kamchatka  peninsula  [see  Kausel 
et  al.  (1977)  for  a discussion  of  this  event].  Analysis  of  this  earthquake  indicates  ( 
that  the  source  is  approximated  by  purely  dip-slip  motion  on  a vertical  fault  plane.  ^ 
We  used  the  station  at  HNR,  which  provides  a N-S  path,  so  that  the  recording  on  \ 
the  E-W  instrument  can  be  considered  to  contain  only  the  azimuthal  component 
of  displacement.  The  focal  depth  is  in  the  intermediate  range.  JThe'initial  earth 
model  that  we  used  in  our  calculations  was  a somewhat  smoothed  version  of  the 
average  oceanic  structure  described  in  Kausel  et  a/.  (1977),  which  is  essentially  the 
CIT-11  structure  (Figure  10,  Toksoz  and  Anderson,  1966).  About  200  layers  were 
used  to  model  the  vertical  heterogeneity  of  the  crust-mantle  system;  an  adaptation 
of  model  MM8  (Anderson  et  al.  1965)  was  used  to  approximate  the  dependence  of 
the  intrinsic  attei^ation  on  depth. 

Tfc  prepare' the  "experimental  and  theoretical  seismograms  for  comparison,  the 
usum  ramp  function  was  removed  from  the  experimental  time  series  (James  and 
Linde,  1971).  The  result  of  this  is  the  upper  trace  in  Figure  1.  Since  the  theoretical 
seismograms  do  not  contain  periods  below  10  sec,  we  ran  a low-pass  frequency 
(boxcar)  filter,  with  cutoff  at  10  sec,  over  the  experimental  data  (second  trace  in 
Figure  1).  A Gaussian  roll-off  filter  with  peak  at  15  sec  and  a 90  per  cent  decay  at 
10  sec,  was  then  applied  to  suppress  the  time-domain  oscillations  introduced  by 
the  abrupt  cutoff  of  the  boxcar  filter.  The  result  of  these  operations  on  the 
experimental  seismogram  is  shown  in  the  third  trace  of  Figure  1.  The  same  Gaussian 
roll-off  filter  is  applied  to  the  theoretical  seismograms,  so  that  both  the  theoretical 
and  experimental  time  series  are  processed  in  the  same  way.  Thus,  we  can  compare 
these  time  series  directly. 

In  our  sample  computations,  we  have  varied  the  source  parameters  over  the 
possible  range  of  variation  from  our  point-source  .solution.  The  first  figure  shows 
the  dip  angle  dependence.  With  the  theoretical  seismograms  available  for  compari- 
son, we  can  easily  identify  several  of  the  arrivals  on  the  experimental  record.  The 
most  striking  of  these  are  the  body-wave  phases  S,  sS.  ScS,  sScS,  and  the  later, 
surface-wave  arrivals.  (The  various  arrivals  are  identified  in  Figure  5.)  For  this 
particular  event,  the  surface  waves  are  composed  of  multimode,  oceanic  Sa  (Kausel 
et  al.,  1977).  In  Figure  1 we  have  varied  the  dip  angle  from  that  of  a vertical  (90°) 
fault  plane,  to  70°;  the  results  show  that  as  dip  angle  decrea.ses,  the  normalized  S 
and  ScS  phases  remain  unchanged,  but  sS  and  sSc5  diminish.  Thus,  the  dip  angle 
ap()ears  to  have  its  main  effect  on  body  waves  when  they  leave  the  focus  above 
the  equator  of  the  focal  sphere.  Both  the  long-  and  short-period  energy  increase,  in 
the  surface-wave  portion  of  the  seismogram,  when  the  dip  angle  decreases  from 
vertical.  Figure  2 illustrates  the  dependence  of  the  seismogram  on  source  depth. 
Here  we  have  varied  this  depth  from  220  to  140  km;  the  consequent  sep.aration  of 
body-wave  arrivals,  and  development  of  the  Sa  wave  train  in  the  lower  set  of 
theoretical  traces,  clearly  exhibit  the  power  of  complete  theoretical  .seismograms 
as  interpretive  tools.  As  Figure  3 illustrates,  the  theoretical  seismogiaras  from  this 
class  of  sources  are  less  sensitive  to  changes 'in  slip  angle  than  to  changes  in  fault 
dip  or  source  depth. 

The  azimuthal-component  seismogram  is  composed,  in  principle,  of  two  contri- 
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butions:  that  from  the  torsi(^l.m^qde_Mciution,  and  that  Jrom.  the  spheroidal' , (.  . ' , 

modes.  The  general  depressions  for  the  polar  (0)  and  azimuthal  («)  components  of.'  J 
displacement  are  given  below  for  the  torsional  (D  and  spheroidal  (S)  excitations  ' 
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Fig.  1.  D^pendenc^  of  computed  sel^moKrams  on  dip  angle  it  llio  vmi  •>  For  cinirast.  the  eaperi-  i 
mentnl  traces  extend  further  to  the  left  and  nghi,  anil  the  diimmant  leastwordi  lube  of  ihe  first  j 
arrival — S — is  upward:  on  the  shorter,  theoretjc.al  seLsmograms,  ea.stward  is  down.  Kach  tra<e,iieTiV)fmaTW 
iied  relative  to  it-s  maximum  di.splacemeni  from  zero.  From  the  top  down,  the  experinienlaFtracrs  slvw  V 
the  E-W  long-penod  seLsmngram  at  HNH  with  the  ramp  function  removed,  then  with  a low. pass  ^ 
frequency  Filter  applied  (cutoff  at  0.1  Hzl  after  ramp  removal,  and  the  lowerawo-aUw-hKvi*  a Gaussian 
roll-off  Filter  applied.  In  Figures  I to  1 we  show  the  effects  of  varying  .source  parameters  from  an  initial 
,.,»..hs..i.  m;  'ip  angle  of  the  faiill  plane  holow  th-  hniaront.al.  <.  of  Qi)'’  azimuthal  angle  between 
strikejjne  and  epiceiiler-to-.siution  line,  of  40';  s-nirce  depth,  h,  Ls  initially  set  .at  IHO  km  iK.ausel  et 
al.  1977,  give  the  reason  fur  this  initial  choice  of  source  depthi;  slip  angle  of  the  hanging  wall  relative  lo 
the  foot  wall,  .\.  of  90°;  and.  a step-hinction  displacement  dislocation  is  a.ssiimed  at  the  focus,  where  a 
point  source  Is  as.siimed.  as  well  as  a component  of  stress,  normal  to  Ihe  fault  plane,  which  is  continiioiis 
dunng  the  dislocation.  [See  Figure  2.  Ksosel  et  of.  1 19771  hu  fault  plane  geomi'irv  ) In  the  al'"ve  rigure. 
the  a.ssocmted  dip  angle  is  given  at  the  left  of  each  theoretical  trace;  the  upper  c.r.iup  of  'hesc  tiiov 
series  represents  actual  grnunfl  displacement,  and  the  lower  set  includes  insinimental  response  [the  10  / 

100  instrument  descnbed  by  KnuiiofTef  of.  1197311  and  Gaussian  ruil-uff  Filter.  . / 
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with  m S 2.  Since  ( u^)t  and  ( u»)s  are  of  the  same  order  of  magnitude,  so  must  V'r 
and  V[s.  Thus,  in  terms  of  an  order-of-magnitude  estimate,  the  relative  sizes  of  the 
two  contributors  to  azimuthal-component  seismograms,  lu^ir  and  (u,>).s-,  are  gov- 
erned by  the  relative  ff  dependences  of  these  contribution:;.  Since  these  dependences 
for  (Ua).s  and  {u^W  are  the  same,  the  order-of-magnitude  comparison  of  the  two 
contributors  to  the  azimuthal  displacement  is  the  same  as  the  comparison  of  the 
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Fig.  2.  Dependence  of  computed  s«-i.smoBriuns  on  source  depth,  which  i.s  given  at  the  left  of  each 
theoretical  setsmogram.  See  caption  of  Figure  I tor  further  deiaila. 


polar  and  azimuthal  di.  plact'inents  torsional  excitation,  i.e. 
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The  results  of  our  direct  computation  of  (u«)r  and  (Uj)r  are  .shown  in  Figtire  d. 
where  we  see  that  relative  to  the  azimuthal  displacements,  the  polar  .seismogram  is 
a straight  line.  In  fact,  the  numerical  results  show  that  the  ratio  of  maximum 
displacements  in  the  two  ca-ses  is  about  lOl  to  1.  We  therefore  conclude  that  for 
normal  WWSSN  seismograms  which  contain  all  of  the  generated  energy  down  to  a 
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FlO.  3.  Dependence  of  compute*)  "eismoiiTam*  *>n  ^lip  angle,  which  i.i  given  n the  left  of  each  theoretical 
seiatnogTam.  See  caption  of  Figure  I for  further  details. 
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periiid  of  10  sec,  the  total  azimuthal  component  of  displacentent  can  be  computed 
by  considering  only  the  torsional-wave  contnbution. 

Before  going  into  a specific,  direct  compjiri.son  of  theoretical  tind  

aeiamugrama,  it  should  be  emphasized  thit  we  have  nut  made  any  .special  aiteinpi 
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to  fit  the  experimental  data  yet.  We  have  just  picked  the  best  fit  to  the  experimental 
record  from  among  the  theoretical  traces  shown  in  Figures  1 to  3.  Even  so,  the 
feature-by- feature  agreement  is  striking  over  much  of  the  record.  This  is  shown  in 
Figure  5.  The  first  point  to  note  is  the  success  of  the  theoretical  work  in  generating 
the  phases  5cS  and  sSc5  that  are  reflected  from  the  core-mantle  boundary.  This 
means  that  we  can  expect  this  kind  of  analysis  to  be  useful  over  the  entire  mantle 
down  to  the  core.  The  second  point  is  the  phese-bv-phase  agreement  in  the  boHv 
waves  identified  in  the  figure,  and  the  de«  r»*ase  in  the  quality  of  the  agu.*emeiii 
lieyond  them.  According  to  the  recent  paper  Ijy  K.au.sel  e.t  nl.  (19771,  the  properties 
of  Sa  are  determined  by  the  structure  from  the  vicinity  of  the  “650-km"  discontinuity 
upward.  Therefore,  the  good  agreement  in  initial  body  waves,  and  the  relatively 
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Fic..  6.  Test  of  relaiiv»  effects  on  initial  body-wave  ph.-i-es  and  on  Inter.  <h  illower  namnlioi:  nrnv.nls 
composing  the  theoretical  seismogram,  when  the  structure  .ihove  a depth  of  rig)  km  is  varied,  i'pper 
and  lower  traces  reproduce  the  experimental  record;  shorter  traces  are  theoretical  sei.smograms  from 
indicated  structure.  These  velocity  structure.s  are  variation.s  from  CIT-11  iKingwoud.  ISTo);  all  u.se  the 
attenuation  model  for  the  Pacific  that  la  given  by  Mitchel  UU76). 

poorer  agreement  in  surface  waves  implies  Mint  we  have  essentially  the  right 
structure  below  7lH)  km,  hut  that  we  need  .sonie  nii"!ilii  oiion  our  eartfi  model 
above  tliat.  As  iiieniioiied  .diove.  the  .structure  we  used  is  an  average  oceanic  model;' 
however,  the  Kamchatka- KNR  path  is  along  the  “old”  edge  of  the  Pacific  plate. 
Recent  seismic  studies  (e.g.,  Leed.s,  et  al.  1974;  Leeds.  107.'))  indicate  that  the 
thickness  of  the  lid  grows  with  increasing  distance  from  the  ridges.  Tlierefore. 
along  the  paths  the  lid  should  be  thicker  than  average.  This  may  be  one  .source  of 
the  discrepancies  between  the  theoretical  and  experimental  traces  in  Figure  5. 
Another  may  be  oceanic  crustal  thicknesses,  for  perhaps  lOtX)  km  north  of  HNR, 
that  may  exceed  30  km  (Furumoto  et  al.,  1973). 

To  give  some  indication  of  whether  modifications  in  the  structure  above  700  km 
can  be  expected  to  improve  the  comparison  in  Figure  5,  we  have  repeated  the 
computations  with  three  modified  structures.  The  structures  and  results  are  show.i 
in  Figure  6,  which  indicates  that  quite  a bit  can  be  done  with  the  single-structure 
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interpretation  before  making  the  obvious  generalization  of  the  technique  to  include  __ 

lateral  heterogeneity.  Consideration  of  the  e.\perimental  coda  is  beyond  the  scope 
of  thisjetter^^  . ^ 


Conclusions 


It  i.s  now  feasible  to  make  direct  comparisons  of  theoretical  and  experimental 
seisni  >grj>fns,  for  toi.sionnl-wave  records  obtained  from  the  long-peiiod  instruments 
of  the  WWSSN.  With  our  |)resent  technique,  all  of  the  generated  energy — at  periods 
above  10  .sec — is  exhibited  on  the  calculated  lime  series,  body  waves  as  well  as 
surface  waves.  Applifation  of  this  result  might  be  expected  to  provide  the  highest 
source  and  structural  resolution  that  has  yet  been  achieved. 

Our  latest  work  (Schwab  el  al.,  1977)  with  spheroidal  waves  on  spherical,  gravi- 
tating models  of  the  earth,  shows  that  even  when  we  employ  highly  opiimired 
versions  of  the  current  algorithms,  these  calculations  are  sLx  times  slower  than  the 
comparable  torsional-wave  computation!!.  We  conclude  from  this  that  new  aleo- 
rithms  are  required  before  the  efficiency^  of  spheroidal-wave  computations  (includ- 
ing gravity)  will  approach  a level  justifying  computation  of  the  associated  theoretical 
seismograms,  down  to  a period  as  low  as  10  sec. 
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